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Introduction 


Background 

An  important  component  of  most  coastal  and  ocean  engineering  projects  is 
an  accurate  assessment  of  the  wave  climate  at  the  project  site.  Typical  applica¬ 
tions  include  determination  of  siltation  rates  inside  entrance  channels  and  harbor 
basins,  determination  of  safe  conditions  for  the  loading/offloading  of  ships,  opti¬ 
mization  of  harbor  layouts  for  both  wind-generated  and  long-period  infragravity 
waves,  design  of  structures  such  as  breakwaters,  and  the  evaluation  of  the  impact 
of  coastal  structures  on  adjacent  shorelines.  Nearshore  wave  conditions  are  norm¬ 
ally  determined  from  deepwater  conditions  because  long-term  wave  data  are 
usually  unavailable  for  most  project  sites.  These  offshore  wave  characteristics 
have  to  be  transformed  to  the  project  site  taking  into  account  the  effects  of  wind- 
wave  generation,  shoaling/refraction  over  seabed  topography,  energy  dissipation 
due  to  wave  breaking  and  bottom  friction,  wave  reflection/diffraction  near  struc¬ 
tures,  nonlinear  wave-wave  interactions,  and  wave  interaction  with  current  fields. 

A  number  of  mathematical  models  have  been  developed  to  simulate  the  prop¬ 
agation  and  transformation  of  waves  in  coastal  regions  and  harbors.  The  different  ‘ 
models  are  based  on  different  assumptions,  which  limit  the  types  of  problems  to 
which  they  can  be  applied.  Examples  include  spectral  wind-wave  models  for 
wave  propagation  in  open  water  where  the  processes  of  wind  input,  shoaling,  and 
refraction  are  dominant;  parabolic  mild-slope  equation  models  for  wave  propaga¬ 
tion  over  large  coastal  areas  when  reflection  is  negligible;  Helmholtz  equation 
models  for  wave  agitation  and  harbor  resonance  in  water  of  constant  depth;  ellip¬ 
tic  mild-slope  models  for  wave  agitation  and  harbor  resonance  in  water  of  vary¬ 
ing  depth;  and  Boussinesq  models  for  nonlinear  wave  refraction-diffraction  in 
shallow  water. 

Numerical  models  available  at  the  U.S.  Army  Corps  of  Engineers  for  pre¬ 
dicting  wave  conditions  in  coastal  regions  and  harbors  include  the  spectral  wind- 
wave  model  STWAVE  (Smith,  Sherlock,  and  Resio  2001)  and  the  elliptic  mild- 
slope  model  CGWAVE  (Demirbilek  and  Panchang  1998).  STWAVE  is  a  wind- 
wave  propagation  model  based  on  the  wave  action  conservation  equation.  It  is  a 
phase-averaged  model,  i.e.,  it  assumes  that  phase-averaged  wave  properties  vary 
slowly  over  distances  of  the  order  of  a  wavelength.  This  allows  the  efficient  com¬ 
putation  of  wave  propagation  over  large  open  coastal  areas.  Due  to  the  phase¬ 
averaging  procedure,  STWAVE  cannot  accurately  resolve  rapid  variations  that 
occur  at  subwavelength  scales  due  to  wave  reflection/diffraction. 
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Phase  resolving  models  based  on  either  the  mild-slope  equation  or 
Boussinesq-type  equations  are  better  suited  for  problems  involving  the 
reflection/difffaction  of  waves  such  as  in  coastal  entrances  and  harbors.  The 
mild-slope  and  Boussinesq  equations  are  vertically  integrated  equations  for  wave 
propagation  in  the  two-dimensional  horizontal  plane  with  different  assumptions 
made  for  the  variation  of  fluid  motion  over  the  water  depth.  The  mild-slope 
equation  derivation  assumes  a  hyperbolic  cosine  variation  of  the  velocity  poten¬ 
tial  over  depth,  consistent  with  linear  monochromatic  waves  in  water  of  arbitrary 
depth,  while  the  Boussinesq  equation  derivation  assumes  a  quadratic  profile, 
valid  for  shallow-water  waves  with  wavelengths  much  longer  than  the  water 
depth. 

This  report  describes  BOUSS-2D,  a  comprehensive  numerical  model  based 
on  a  time-domain  solution  of  Boussinesq-type  equations.  The  classical  form  of 
the  Boussinesq  equations  for  wave  propagation  over  water  of  variable  depth  was 
derived  by  Peregrine  (1967).  The  equations  were  restricted  to  relatively  shallow- 
water  depths,  i.e.,  the  water  depth,  h,  had  to  be  less  than  one-fifth  of  the  wave¬ 
length,  L,  in  order  to  keep  errors  in  the  phase  velocity  to  less  than  5  percent. 
Nwogu  (1993)  extended  the  range  of  applicability  of  Boussinesq-type  equations 
to  deeper  water  by  recasting  the  equations  in  terms  of  the  velocity  at  an  arbitrary 
distance  za  from  the  still-water  level,  instead  of  the  depth-averaged  velocity.  The 
elevation  of  the  velocity  variable  za  becomes  a  free  parameter,  which  is  chosen  to 
optimize  the  linear  dispersion  characteristics  of  the  equations.  The  optimized 
form  of  the  equations  has  errors  of  less  than  2  percent  in  the  phase  velocity  from 
shallow-water  depths  up  to  the  deepwater  limit  (h/L  =  0.5). 

Despite  the  improvement  in  the  frequency  dispersion  characteristics, 
Nwogu’s  (1993)  equations  are  based  on  the  assumption  that  the  wave  heights 
were  much  smaller  than  the  water  depth.  This  limits  the  ability  of  the  equations 
to  describe  highly  nonlinear  waves  in  shallow  water  and  led  Wei  et  al.  (1995)  to 
derive  a  fully  nonlinear  form  of  the  equations.  The  fully  nonlinear  equations  are 
particularly  useful  for  simulating  highly  asymmetric  waves  in  shallow  water, 
wave-induced  currents,  wave  setup  close  to  the  shoreline,  and  wave-current 
interaction. 

As  ocean  waves  approach  the  shoreline,  they  steepen  and  ultimately  break. 
The  turbulence  and  currents  generated  by  breaking  waves  are  important  driving 
mechanisms  for  the  transport  of  sediments  and  pollutants.  Nwogu  (1996) 
extended  the  fully  nonlinear  form  of  the  Boussinesq  equations  to  the  surf  zone, 
by  coupling  the  mass  and  momentum  equations  with  a  one-equation  model  for 
the  temporal  and  spatial  evolution  of  the  turbulent  kinetic  energy  produced  by 
wave  breaking.  The  equations  have  also  been  modified  to  include  the  effects  of 
bottom  friction  and  flow  through  porous  structures.  The  modified  equations  can 
simulate  most  of  the  hydrodynamic  phenomena  of  interest  in  coastal  regions  and 
harbor  basins  including: 

a.  Shoaling. 

b.  Refraction. 

c.  Diffraction. 
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d.  Full/partial  reflection  and  transmission. 

e.  Bottom  friction. 

/  Nonlinear  wave-wave  interactions. 

g.  Wave  breaking  and  runup. 

h.  Wave-induced  currents. 

i.  Wave-current  interaction. 

A  time  domain,  finite-difference  method  is  used  to  solve  the  Boussinesq 
equations.  The  area  of  interest  is  discretized  as  a  rectangular  grid  with  the  water- 
surface  elevation  and  horizontal  velocities  defined  at  the  grid  nodes  in  a  stag¬ 
gered  manner.  Time-histories  of  the  velocities  and  fluxes  corresponding  to  inci¬ 
dent  storm  conditions  are  specified  along  external  or  internal  wave  generation 
boundaries.  The  wave  conditions  may  be  periodic  or  nonperiodic,  unidirectional 
or  multidirectional.  Waves  propagating  out  of  the  computational  domain  are 
absorbed  in  damping  layers  placed  around  the  perimeter  of  the  numerical  basin. 
Porosity  layers  can  be  placed  inside  the  computational  domain  to  simulate  the 
reflection  and  transmission  characteristics  of  structures  such  as  breakwaters. 


Purpose 

This  report  describes  BOUSS-2D,  a  comprehensive  numerical  model  for 
simulating  the  propagation  and  transformation  of  waves  in  coastal  regions  and 
harbors  based  on  a  time-domain  solution  of  Boussinesq-type  equations.  An 
overview  of  the  theoretical  background  behind  the  model  is  described  in 
Chapter  2.  The  numerical  scheme  used  to  solve  the  equations  is  described  in 
Chapter  3.  The  steps  involved  in  setting  up  and  running  the  model  are  described 
in  Chapter  4.  A  number  of  analytical,  laboratory,  and  field  test  cases  have  also 
been  used  to  validate  the  model  in  Chapter  5.  The  different  test  cases  were 
selected  in  Chapter  5  to  evaluate  the  ability  of  the  model  to  deal  with  individual 
wave  transformation  processes  such  as  refraction,  diffraction,  wave  breaking, 
etc.,  as  well  as  the  combination  of  processes  that  occur  in  practical  engineering 
problems.  Chapter  6  provides  a  summary  and  conclusions. 
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2  Theoretical  Background 


Governing  Equations 

B0USS-2D  is  based  on  Boussinesq-type  equations  derived  by  Nwogu  (1993, 
1996).  The  equations  are  depth-integrated  equations  for  the  conservation  of  mass 
and  momentum  for  nonlinear  waves  propagating  in  shallow  and  intermediate 
water  depths.  They  can  be  considered  to  be  a  perturbation  from  the  shallow-water 
equations,  which  are  often  used  to  simulate  tidal  flows  in  coastal  regions.  For 
short-period  waves,  the  horizontal  velocities  are  no  longer  uniform  over  depth 
and  the  pressure  is  nonhydrostatic.  The  vertical  profile  of  the  flow  field  is 
obtained  by  expanding  the  velocity  potential,  O,  as  a  Taylor  series  about  an 
arbitrary  elevation,  za,  in  the  water  column.  For  waves  with  length,  L,  much 
longer  than  the  water  depth,  h,  the  series  is  truncated  at  second  order  resulting  in 
a  quadratic  variation  of  the  velocity  potential  over  depth: 

<D(x,z,0  =  <j»„  +\i2{za -z)[V<)>a  -VA] 

1,2  (1) 

+  ^-[(za  +  h)2  -  (z  +  h)2 ] v2(|)0  +  0(p4) 

where  (f>a  =  0(;c,za,?),  V  =  (9/9 x,d/dy),  and  (J.  =  h/L  is  a  measure  of  frequency 
dispersion.  The  horizontal  and  vertical  velocities  are  obtained  from  the  velocity 
potential  as: 

«a  +  (Za  -Z)[V(«a  ' V K )  +  (V  *  «a)V^] 

i[(z„  +  *)!-(z  +  /1)!]v(V.«J  (2) 

-[ua-Vh  +(z  +  h)V-ua]  (3) 


u(x,z,t )  =  VO  = 


+ 


,  ,  90 

w(x,z,t)  =  —  = 


where  ua  =  VO|Za  is  the  horizontal  velocity  at  z  =  za.  Given  a  vertical  profile  for 
the  flow  field,  the  continuity  and  Euler  (momentum)  equations  can  be  integrated 
over  depth,  reducing  the  three-dimensional  problem  to  two  dimensions.  For 
weakly  nonlinear  waves  with  height,  H,  much  smaller  than  the  water  depth,  h,  the 
vertically  integrated  equations  are  written  in  terms  of  the  water-surface  elevation 
T \(x,t)  and  velocity  ua(x,t)  as  (Nwogu  1993): 


4 


Chapter  2  Theoretical  Background 


0 


(4) 


T)(  +  V  •  Ify  = 

«a ,/ + gVn + («a  •  v)«a  +  za  [v  («a ,  •  v/?) + (V  ■  «a>,)V/z] 


(5) 


where  g  is  the  gravitational  acceleration  and  «/is  the  volume  flux  density  given 
by: 


«/  =  \\“dz  =  (A +Tfl  «„  +  A^z0  + 1  j[V(«a  •  Vh) + (V  •  kJVA] 


+ 


(^a  +  A)2 


(6) 


V(V-«a) 


The  depth-integrated  equations  are  able  to  describe  the  propagation  and  trans¬ 
formation  of  irregular  multidirectional  waves  over  water  of  variable  depth.  The 
elevation  of  the  velocity  variable  za  is  a  free  parameter  and  is  chosen  to  minimize 
the  differences  between  the  linear  dispersion  characteristics  of  the  model  and  the 
exact  dispersion  relation  for  small  amplitude  waves.  The  optimal  value, 
za  =  -0.535A,  is  close  to  middepth. 

For  steep  near-breaking  waves  in  shallow  water,  the  wave  height  becomes  of 
the  order  of  the  water  depth  and  the  weakly  nonlinear  assumption  made  in  deriv¬ 
ing  Equations  4  and  5  is  no  longer  valid.  Wei  et  al.  (1995)  derived  a  fully  non¬ 
linear  form  of  the  equations  from  the  dynamic  free  surface  boundary  condition  by 
retaining  all  nonlinear  terms,  up  to  the  order  of  truncation  of  the  dispersive  terms. 
Nwogu  (1996)  derived  a  more  compact  form  of  the  equations  by  expressing  some 
of  the  nonlinear  terms  as  a  function  of  the  velocity  at  the  free  surface,  u^,  instead 
of  ua.  Additional  changes  have  also  been  made  to  the  equations  to  allow  for 
weakly  rotational  flows  in  the  horizontal  plane  and  ensure  that  za  remains  in  the 
water  column  for  steep  waves  near  the  shoreline  and  during  the  wave  runup 
process.  The  revised  form  of  the  fully  nonlinear  equations  can  be  written  as: 

Tl+V-n,  =  0  (7) 

«<M  +  &Vrl  +  (»„ -VK  +  +  (z0  -T|)[v  («„,  •  VA)  +  (V  •  «„  ,)VA] 

-  [(«lM  VA)+(A+TDV.«tM]vn 

+  [V(**0>,  •  VA)  +  (V-«a,)VA  +  (za  +  A)V (V-«a)]za(  =  0 


Chapter  2  Theoretical  Background 


5 


where  za  is  now  a  function  of  time  and  is  given  by  za+h  -  0.465(/i+'n).  The 
volume  flux  density  «/  is  given  by: 


«/  =  (^+1l)i««  + 


(z„+A)- 


(A  +  ri) 


[V(«„  •VA)+(V-«0)VA] 


(za+h)2  (h  +  Ti)2 


V(V«a) 


(9) 


The  fully  nonlinear  equations  are  able  to  implicitly  model  the  effects  of  wave- 
current  interaction.  Currents  can  either  be  introduced  through  the  boundaries  or 
by  explicitly  specifying  a  current  field,  U. 


Linear  Dispersion  Properties 

The  linear  dispersion  relation  of  the  Boussinesq  model  that  relates  the 
wavelength,  L,  to  the  wave  period,  T,  is  given  by  Nwogu  (1993)  as: 


l-(q  +  l/3 Xkhf 
T2~[  1  -  a(kh)2 


(10) 


where  C  is  the  phase  speed,  k  =  2 n/L  is  the  wave  number,  and  a  = 

[( za+h)2/h 2  - 1]/2.  Depending  on  the  elevation  of  the  velocity  variable  or  the  value 
of  a,  different  dispersion  relations  are  obtained.  If  the  velocity  at  the  seabed 
( za  =  -h)  is  used,  a  =  -1/2.  Alternatively,  if  the  velocity  at  the  still-water  level 
(za  =  0)  is  used,  a  =  0.  The  dispersion  relation  of  the  classical  form  of  the 
Boussinesq  equations  which  uses  the  depth-averaged  velocity  as  the  velocity 
variable  corresponds  to  a  =  -1/3.  Witting  (1984)  obtained  the  value  a  =  -2/5  from 
the  Pade  (2,2)  approximant  of  tanh  kh. 

The  phase  speeds  for  different  values  of  a,  normalized  with  respect  to  the 
linear  theory  phase  speed  are  plotted  as  a  function  of  relative  depth  in  Figure  1 . 
The  relative  depth  is  defined  as  the  ratio  of  the  water  depth,  h,  to  the  equivalent 
deepwater  wavelength  L0  -  gT2 12 k.  The  deepwater  depth  limit  corresponds  to 
h/L  =  0.5.  The  different  dispersion  equations  are  all  equivalent  in  relatively 
shallow  water  ( h/L  <  0.02/  but  gradually  diverge  from  the  exact  solution  with 
increasing  depth.  An  optimal  depth  za  =  -0.535 h  gives  errors  of  less  than  2  per¬ 
cent  in  the  phase  velocity  from  shallow-water  depths  up  to  the  deepwater  limit. 


Nonlinear  Properties 

In  an  irregular  sea  state,  different  frequency  components  interact  to  generate 
forced  waves  at  the  sum  and  difference  frequencies  of  the  primary  waves  because 
of  the  nonlinear  nature  of  the  boundary  condition  at  the  free  surface.  Consider  a 


6 


Chapter  2  Theoretical  Background 


(1994)  derived  similar  expressions  for  the  weakly  nonlinear  form  of  the 
Boussinesq  equations: 


G±((apco2) 


{0,cq2  {Khf  cos  A9[l  -  (0  +  1/ 3)(k±hf  ] 

2  Xk[k'2h3 

co±[l  -  a.(k±h)2][(i)1k2h(klh  ±  k2h  cos  AQ)  +  co2At1/i(A:1/i  cosAQ  ±  k2hj\ 

2  Xk[k'2h3 


(13) 


where Dq  =  q\  ±  q2,  k ±  =  |  k\±  k2\y  k'  =  k  [1  -  (a +  1/3)  (kh) 2],  and 
A,  =  eo^[l  -  a (k±h)2 ]  -  gk2±h[l  -  (a  + 1  / 3 )(k±h)2  ] 


(14) 


Figure  2  shows  a  comparison  of  the  quadratic  transfer  function  of  the  weakly 
nonlinear  Boussinesq  model  with  that  of  second-order  Stokes  theory  for  unidirec¬ 
tional  waves  where  the  wave  group  period  is  10  times  the  average  of  the  indi¬ 
vidual  wave  periods,  i.e.,  co_  =  (tOj  +  o)2)/20.  The  weakly  nonlinear  Boussinesq 
model  underestimates  the  magnitude  of  the  setdown  wave  and  second  harmonic 
at  the  deepwater  depth  limit  by  65  percent  and  45  percent  respectively.  Hence,  it 
cannot  accurately  simulate  nonlinear  effects  in  deep  water.  To  reasonably  simu¬ 
late  nonlinear  effects,  the  weakly  nonlinear  model  should  be  restricted  to  the 
range  0  <  h/L  <  0.3. 


Simulation  of  Wave  Breaking 

The  turbulent  and  highly  rotational  flow  field  under  breaking  waves  is 
extremely  complex  and  difficult  to  model  even  with  the  Reynolds-averaged  form 
of  the  Navier-Stokes  equations  (e.g.,  Lin  and  Liu  1998;  Bradford  2000).  In 
BOUSS-2D,  we  do  not  attempt  to  model  details  of  the  turbulent  motion,  but 
rather,  simulate  the  effect  of  breaking-induced  turbulence  on  the  flow  field.  We 
have  tried  to  develop  a  generic  model  that  can  be  applied  to  regular  or  irregular 
waves,  unidirectional  or  multidirectional  waves,  and  simple  or  complex  bottom 
topography  without  having  to  recalibrate  the  model  each  time.  The  key  assump¬ 
tions  made  in  developing  the  model  are  (see  Nwogu  1996): 

a.  The  breaking  process  is  assumed  to  be  “spilling.” 

b.  Turbulence  is  produced  in  the  near-surface  region  when  the  horizontal 
velocity  at  the  free  surface,  «n,  exceeds  the  phase  velocity,  C. 

c.  The  rate  of  production  of  turbulent  kinetic  energy  is  proportional  to  the 
vertical  gradient  of  the  horizontal  velocity  at  the  free  surface,  9«/3z|z=n. 

d.  Breaking-induced  turbulence  is  primarily  convected  in  the  near-surface 
region  with  the  horizontal  velocity  at  the  free  surface. 
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wave  train  consisting  of  two  small  amplitude  periodic  waves  with  amplitudes,  a\ 
and  a2,  frequencies,  coj  and  ©2,  wave  numbers,  k{  and  k2,  and  propagating  in 
directions  0]  and  02  respectively.  The  water-surface  elevation  can  be  written  as: 

T)(1)(jc,t)  =  a,  cos(kl  ■  x  -  C0j?)  +  a2  cos (k2  •  x  -  (02t)  (1 1) 


where  k  =  (kcosO,  foin0).  The  second-order  wave  will  consist  of  a  subharmonic 
at  the  difference  frequency,  ©_=©;-  co2,  and  higher  harmonics  at  the  sum  fre¬ 
quencies  2©j,  2cp2  and  ©+  =  ©1  +  ©2.  This  can  be  written  as: 

rf>  (x,  t)  =  -±-  G+  (©j ,  ©! ,  0( ,  0j )  (cos(2^!  •  x  -  2©jt) 

2 

+  — -G+(©2,©2,02>02)  cos(2k2- x-2(o2t)  (12) 

+  cos(2*±  •  x - 2©±t) 


where  k+  =  k\±  k2,  and  G±(©],  ©2,  0i,02)  is  a  bidirectional  quadratic  transfer  func¬ 
tion  that  relates  the  amplitude  of  the  second-order  waves  to  the  amplitude  of  the 
first-order  waves.  Dean  and  Sharma  (1981)  derived  expressions  for  the  bidirec¬ 
tional  quadratic  transfer  function  based  on  second-order  Stokes  theory.  Nwogu 
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Figure  2.  Comparison  of  quadratic  transfer  function  for  Boussinesq  and  Stokes 
theories 

The  effect  of  wave  energy  dissipation  due  to  breaking  is  simulated  in  the 
Boussinesq  model  by  introducing  an  eddy  viscosity  term  to  the  right-hand  side  of 
the  momentum  equation  (Equation  5  or  8).  Nwogu  (1996)  used  a  dissipative  term 
of  the  following  form: 

KeMng  =  -v,V(V.«J  (15) 

where  v,  is  the  turbulent  eddy  viscosity.  As  pointed  out  by  Kennedy  et  al.  (2000), 
it  is  important  for  the  dissipative  term  to  dissipate  energy  but  conserve  momen¬ 
tum  to  accurately  capture  details  of  the  mean  flow  field  associated  with  breaking 
waves.  A  modified  form  of  the  dissipative  term  that  ensures  that  momentum  is 
conserved  can  be  written  as: 


1 

A+rj 


V{v,(/z+Ti)V-Ma} 


(16) 


The  eddy  viscosity  is  determined  from  the  amount  of  turbulent  kinetic  energy,  k, 
produced  by  wave  breaking,  and  a  turbulence  length  scale,  using: 


V,  =  yfkl, 


(17) 
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A  one-equation  model  is  used  to  describe  the  production,  advection,  diffusion, 
and  dissipation  of  the  turbulent  kinetic  energy  produced  by  wave  breaking: 


k.  — 


■u^Vk  +  oV  •  V  (v,£)  +  B 


Jc~o 


2 

1 

2  “ 

.w 

t 

KdZj 

2=11 


-  C„—  (18) 


The  waves  are  assumed  to  start  breaking  when  the  horizontal  component  of  the 
orbital  velocity  at  the  free  surface,  un,  exceeds  the  phase  velocity  of  the  waves, 
C.  The  parameter  B  is  introduced  to  ensure  that  production  of  turbulence  occurs 
after  the  waves  break,  i.e., 


B 


o  Ki<c 

i  M^c 


(19) 


The  phase  velocity  is  determined  from  the  linear  dispersion  relation  (Equa¬ 
tion  10)  using  the  average  zero-crossing  period  of  the  incident  wave  train.  While 
this  approach  leads  to  some  waves  in  an  irregular  wave  train  breaking  at  slightly 
different  locations  in  the  model  than  they  would  in  nature,  it  was  found  to  be 
more  stable  than  trying  to  estimate  a  time-dependent  phase  velocity  using 
expressions  such  as  C(t)  =  -r|,/|Vr||. 

The  empirical  constants  CD  and  a  have  been  chosen  as  0.02  and  0.2  respec¬ 
tively.  The  turbulent  length  scale,  /„  remains  the  only  free  parameter  in  the  turbu¬ 
lence  model  and  is  determined  from  comparisons  of  numerical  model  results  with 
experimental  data.  Recommended  values  are  the  significant  wave  height  (/,  = 
Hmo)  for  irregular  waves,  and  the  wave  height  (/,  =  H)  for  regular  waves. 


Bottom  Friction 

The  bottom  boundary  layer  in  wave  fields  is  typically  confined  to  a  tiny 
region  above  the  seabed,  unlike  river  and  tidal  flows  where  it  extends  all  the  way 
up  to  the  free  surface.  There  is,  thus,  very  little  wave  energy  attenuation  due  to 
bottom  friction  over  typical  wave  propagation  distances  of  O  (1km)  used  in 
Boussinesq-type  models.  The  bottom  friction  factor,  however,  plays  a  more 
important  role  in  wave  transformation  close  to  the  shoreline  and  nearshore  circu¬ 
lation  patterns.  The  effect  of  energy  dissipation  due  to  a  turbulent  boundary  layer 
at  the  seabed  has  been  modeled  by  adding  a  bottom  shear  stress  term  to  the  right- 
hand  side  of  the  momentum  equation  (Equation  5  or  8): 


bjriction 


h+T] 


f  u  u 

J  h  a  a 


(20) 


where  fw  is  the  wave  friction  factor.  The  bottom  friction  term  can  also  be  written 
in  terms  of  the  Chezy  coefficient,  Cf,  used  in  tidal  flows  by  replacing/,  with 

g/Cf. 
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Equation  20  has  been  expressed  in  terms  of  «„  instead  of  the  velocity  at  the 
seabed,  ub,  to  minimize  the  additional  computational  expense  of  evaluating  uh. 
The  values  of  the  friction  factors  specified  in  the  model  would  thus  be  slightly 
different  than  those  based  on  the  bottom  velocity. 
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3  Numerical  Solution 


Finite  Difference  Scheme 

The  weakly  and  fully  nonlinear  Boussinesq  equations  (Equations  4-9)  are 
solved  in  the  time-domain  using  a  finite-difference  method.  The  computational 
domain  is  discretized  as  a  rectangular  grid  with  grid  sizes  Ax  and  Ay,  in  the  x  and 
y  directions,  respectively.  The  equation  variables  r|,  ua,  and  va  are  defined  at  the 
grid  points  in  a  staggered  manner  as  shown  in  Figure  3.  The  water  depth  and 
surface  elevation  are  defined  at  grid  points  (ij),  while  the  velocities  are  defined 
half  a  grid  point  on  either  side  of  the  elevation  grid  points.  The  external  boun¬ 
daries  of  the  computational  domain  correspond  to  velocity  grid  points. 


y 

i 


Ajc 


►  X 


Figure  3.  Computational  grid  for  finite  difference  scheme 


12 


Chapter  3  Numerical  Solution 


The  numerical  solution  scheme  is  an  implicit  Crank-Nicolson  scheme  with  a 
predictor-corrector  method  used  to  provide  the  initial  estimate.  The  first  step  in 
the  solution  scheme  is  the  predictor  step  in  which  values  of  the  variables  at  an 
intermediate  time-step  t  =  (n  +  54)At  are  determined  using  known  values  at  t  = 
nAt  (At  is  the  time-step  size).  The  second  step  is  the  corrector  step  in  which  pre¬ 
dicted  values  at  t  =  (n+  54)At  are  used  to  provide  an  initial  estimate  of  the  values 
at  t  =  («+l)At.  The  last  step  is  an  iterative  Crank-Nicolson  scheme,  which  is 
repeated  until  convergence. 

The  partial  derivatives  are  approximated  using  a  forward  difference  scheme 
for  time  and  central  difference  schemes  for  the  spatial  variables.  The  resulting 
Crank-Nicolson  formulation  for  the  weakly  nonlinear  form  of  the  mass  and 
momentum  equations  (Equations  4-6)  can  be  written  as: 


8/n =-  -8®v;+1/2 


(21) 


5, («„  +/A*)  =  -  tftr2 


(«r,2‘)  -^(vr2* 

[s»  +  V(*) 


A 


(22) 


s,  (v0  +  /48<2)VK  +/25^va)  =  -g8'V+1/2  -i84(vr/2')2l-5f  (vr/2‘ 

-  /4[8®W  §,81°  (c177")  +  8®  (h)  8,8^  (23) 

-  fAK<m 

where  the  volume  flux  densities  uf  and  v/are  given  by: 


Ufi+V2,j  =  h  +  T\KiW2J  +  ^/(SxA.+S^) 

+  Ax/3 [28®(A)8f p)  +  8f(h)  8®(^)] 


(24) 


v/u+v2  =  h+r\valJ+U2  +  hy  fx  (8^ua  +  8^va ) 

+  F/3[28<1)(A)8fva+8«W8f(^)  +  8 f(h)  8«(^)] 

The  parameters /i  to  fA  are  given  by: 

f_\(za  +  hf  X 

h  ~  2  6 


(25) 


(26a) 
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(26b) 
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f2  = 


(gg 

2 


/3  = 


(Za.+h)~^ 


(26c) 


/4=[(za +/>)-*] 


(26d) 


The  average  and  difference  operators  in  Equations  21  to  25  are  given  by: 

5,<J>  =  7-(<t>"+1-  f)  (27a) 

A/ x 

f+1/2  =!((()"+  f+1)  (27b) 


<T  =|(<t>,+u+<f>/j)  or  |(t+i/2,y+4>M/2./) 

(27c) 

V  =^(ly+i+(t,u)  or  ^(^1/2+ ^.,-1/2) 

(27d) 

^  =  ^(iu  -40 or 

(27e) 

0r  -bj-mj) 

(27f) 

^  -  2Ax  (<t>/+3/2,Jr  §i-U2,j  ) 

(27g) 

^  ^  ~  2 (^<.2+3/2  ^/.y-i/2 ) 

(27h) 

5.4>  —  ^2  (4Vf-3/2J  2§i+l/2j  +4)/-l/2j) 

(27i) 

5^4)  ”  ^  2  (4)i,y+3/2  24>«J+l/2  +4)|J-l/2) 

(27j) 

(27k) 

Chapter  3  Numerical  Solution 


and  <|>  represents  the  variables  (tj,  ua,  va).  The  finite-difference  operator  defini¬ 
tions  in  Equation  27  corresponding  to  the  x  and  y  components  of  the  momentum 
equation  are  centered  at  (i+l/2j)  and  (i,j+ 1/2),  respectively. 

The  finite  difference  formulation  of  the  continuity  equation  (Equation  21) 
yields  an  algebraic  equation  that  is  explicitly  solved  for  r|  at  all  grid  points.  The 
formulation  for  the  x  andy  momentum  equations  (Equations  22  and  23)  have 
been  decoupled  by  placing  the  vxt,  vyt,  and  v^,  terms  on  the  right-hand  side  of  the 
x  equation  and  the  uxt,  uyh  and  u^,  terms  on  the  right-hand  side  of  the  y  equation. 
This  reduces  the  momentum  equations  to  tridiagonal  equations  for  ua  and  v„ 
along  lines  in  the  x  andy  direction,  respectively.  Tridiagonal  matrices  are  much 
easier  to  store  and  solve  than  the  large  sparse  matrix  equation  that  would  be 
obtained  if  the  equations  were  not  decoupled.  The  major  disadvantage  of  this 
approach,  however,  is  that  the  iterative  step  takes  longer  to  converge  for  shorter- 
period  waves  propagating  at  large  angles  to  the  grid  where  the  higher-order 
cross-derivative  terms  (u^,,  y^,)  become  comparable  in  magnitude  to  the  inline 
derivative  terms  («**,,  v^). 

The  numerical  scheme  is  stable  provided  that  the  Courant  number,  CR,  is  less 
than  1,  i.e., 


Cr 


C2At2 


1  1 

-  +  - 


Ax2  Ay2 


yj 


(28) 


where  C  is  the  phase  speed  based  on  the  average  zero-crossing  period  of  the 
incident  waves.  It  is,  however,  recommended  that  the  Courant  number  be  kept 
within  the  range  0.5  to  0.7  since  nonlinear  wave-wave  interactions,  wave  break¬ 
ing,  and  the  presence  of  reflected  waves  can  affect  the  stability  criterion  of  the 
numerical  model. 


Boundary  Conditions 

To  solve  the  governing  equations,  appropriate  boundary  conditions  have  to 
be  imposed  at  the  boundaries  of  the  computational  domain.  This  requires  specifi¬ 
cation  of  waves  propagating  into  the  domain  and  the  absorption  of  waves  propa¬ 
gating  out  of  the  domain.  The  equations  have  also  been  modified  to  simulate 
wave  interaction  with  fully/partially  reflecting  structures  within  the  computa¬ 
tional  domain.  The  types  of  boundaries  considered  in  BOUSS-2D  include: 

a.  Fully  reflecting  or  solid  wall  boundaries. 

b.  External  wave  generation  boundaries. 

c.  Internal  wave  generation  boundaries. 

d.  Wave  absorption  or  damping  regions. 

e.  Porous  structures. 
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Solid  wall  boundaries 


Along  solid  wall  or  fully  reflecting  boundaries,  the  horizontal  velocity 
normal  to  the  boundary  must  be  zero  over  the  entire  water  depth,  i.e., 

un  =  0  —h  <  z<T|  (29) 


where  n  is  the  normal  vector  to  the  boundary.  This  condition  is  satisfied  in  the 
depth-integrated  equations  by  specifying  that  both  the  volume  flux  density  and 
velocity  normal  to  the  boundary  are  zero,  i.e., 

«„«  =  0  (30) 

ufn  =  0  (31) 

Since  the  equations  are  solved  on  a  staggered  grid,  the  boundary  conditions  are 
specified  as  either  ua  =  Uf  =  0  along  wall  boundaries  perpendicular  to  the  x-axis, 
or  v„  =  v/=  0  along  wall  boundaries  perpendicular  to  the  y-axis. 


External  wave  generation  boundaries 

Along  external  wave  generation  boundaries,  time-histories  of  velocities  ua  or 
va,  and  flux  densities  Uf  or  vf  corresponding  to  an  incident  storm  condition  are 
specified.  The  time-histories  may  correspond  to  regular  or  irregular,  unidirec¬ 
tional  or  multidirectional  waves. 

Regular  Waves.  Regular,  long-crested  wave  conditions  are  specified  in 
terms  of  a  wave  height,  H,  wave  period,  T,  and  direction  of  propagation,  0.  The 
water-surface  elevation  for  small-amplitude  waves  (H«  h )  may  be  written  as: 


H 

—  cos 
2 


(Ax  cos  0  +  ky  sin  0  -  cot) 


(32) 


where  co  =  2n/T  is  the  angular  frequency,  k  is  the  wave  number,  and  0  is  the 
direction  of  wave  propagation  relative  to  the  positive  x-axis.  The  boundary 
conditions  along  a  wave  generation  line  perpendicular  to  the  x-axis  may  be 
obtained  from  the  linearized  form  of  the  continuity  equation  (Equation  5)  for 
water  of  constant  depth  as: 


ua(xg,yg,t)  =  T„(co)cos0  r\(xg,yg,t) 


(33) 


(h  +  T\)  -  h 

ra\ 

m2 

£  o 
^  ) 

_ 

ua(xg>yg’0 


(34) 
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where  (xg,  yg)  are  the  x-y  coordinates  of  the  wave  generation  line  and  Tu( a>)  is  a 
linear  transfer  function  given  by: 


w 


co 


kh 


1- 


^  al 


2  6 


(kh)2 


(35) 


and  a2  =  (za  +h)/h.  Similar  expressions  apply  to  va  and  v/for  wave  generation 
lines  perpendicular  to  the  y-axis  with  cos9  replaced  by  sin0  in  Equation  33. 

For  large  amplitude  waves  [H=  0(h)],  higher  harmonic  components  (2co, 

3co,  etc.)  are  generated  due  to  the  nonlinear  terms  in  the  governing  equations. 
These  waves  are  often  called  bound  waves  since  they  are  attached  to  the  primary 
wave  and  travel  at  its  phase  speed,  C  =  co Ik,  instead  of  the  phase  speed  of  a  free 
wave  at  the  corresponding  frequency.  The  wave  shape  also  changes  from  the 
sinusoidal  shape  assumed  in  Equation  32  to  an  asymmetric  one  with  peaked 
crests  and  broad  shallow  troughs.  If  linear  wave  conditions  are  imposed  at  the 
boundaries,  the  numerical  model  will  generate  free  higher  harmonic  components 
with  the  same  magnitude  but  180  deg  out  of  phase  with  the  bound  waves  at  the 
wavemaker  to  satisfy  the  linear  boundary  condition.  The  presence  of  bound  and 
free  higher  frequency  waves  that  travel  at  different  speeds  will  lead  to  a  spatially 
nonhomogenous  wave  field  with  the  wave  height  and  shape  changing  continu¬ 
ously  over  the  computational  domain. 

The  Fourier  approximation  method  of  Rienecker  and  Fenton  (1981)  has  been 
used  to  solve  the  weakly  nonlinear  form  of  the  Boussinesq  equations  and  develop 
nonlinear  boundary  conditions  for  the  generation  of  large-amplitude  regular 
waves  in  shallow  water.  The  partial  differential  Equations  4  to  6  are  initially 
transformed  into  a  set  of  coupled  nonlinear  ordinary  differential  equations  in 
terms  of  a  moving  coordinate  system,  ^  =  x-Ct.  The  velocity  variable  ua  is 
expanded  as  a  Fourier  series  and  substituted  into  the  governing  equations,  which 
are  evaluated  at  a  finite  number  of  collocation  points  over  half  a  wavelength  to 
yield  a  system  of  nonlinear  algebraic  equations.  A  Newton-Raphson  iterative 
procedure  is  used  to  solve  the  nonlinear  equations  for  the  unknown  values  of  the 
free  surface  displacement  at  the  collocation  points,  the  Fourier  coefficients,  the 
wave  number,  and  the  phase  speed.  Details  of  the  technique  are  provided  in 
Appendix  A. 

Irregular  Unidirectional  Waves.  For  nonperiodic  waves,  the  incident  wave 
conditions  are  typically  expressed  in  the  form  of  a  wave  spectrum,  which 
describes  the  frequency  distribution  of  wave  energy.  Different  parametric  shapes 
have  been  proposed  for  wave  spectra  including: 

a.  The  Pierson-Moskowitz  (PM)  spectrum  for  fully  developed  sea  states  in 
deep  water,  which  is  defined  in  terms  of  the  wind  speed. 

b.  The  Bretschneider  (1959)  spectrum,  which  has  the  same  shape  as  the  PM 
spectrum  but  is  defined  in  terms  of  the  significant  wave  height  and  peak 
period. 
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c.  The  JONSWAP  spectrum  for  fetch-limited  conditions  in  deep  water. 

d.  The  TMA  spectrum  for  fetch-limited  conditions  in  shallow  water. 

e.  The  Ochi-Hubble  spectrum  for  bimodal  sea  states  with  double-peaked 
spectra. 

Expressions  for  the  different  wave  spectra  are  provided  in  Appendix  B. 

The  Fourier  series  technique  is  used  to  generate  time-histories  of  the  velocity 
boundary  conditions  from  the  wave  spectrum.  The  water-surface  elevation,  r|(t), 
is  assumed  to  be  a  zero-mean,  stationary,  random  Gaussian  process.  The  surface- 
elevation  time  series  at  a  reference  point  xr  =  (xr,yr)  in  the  computational  domain 
can  be  represented  as  a  linear  superposition  of  N  regular  wave  components,  i.e., 

N 

n(xr , 0  =  s  ai cos \kJ  ■xr-(Ojt+ ]  (36) 

j= i 

where  ajt  co,-  £,-,  and  kj  =  (kjCOsQ,  kjsinQ)  are  the  amplitude,  angular  frequency, 
phase  angle,  and  wave  number  vector  of  the  fh  frequency  component,  respec¬ 
tively.  The  angle,  0,  is  the  direction  of  wave  propagation  relative  to  the  positive 
x-axis. 

The  wave  spectrum  is  divided  into  N  frequency  bands  with  uniform  spacing, 
Aco,  so  that  the  frequency  of  the  j'h  wave  component  is  given  by  co;  =yAco.  The 
amplitudes  of  the  individual  wave  components  are  obtained  deterministically 
from  the  wave  spectrum,  5n( <o),  as: 

aj  =  ^  (co,)  Aco  (37) 


while  the  phase  angles,  ey,  are  randomly  selected  from  a  uniform  distribution 
between  0  and  2n.  Incident  wave  conditions  are  more  typically  specified  in  terms 
of  the  repeat  period  or  duration  of  the  record,  TD,  and  time-step,  At.  The  values  of 
N  and  Aco  can  be  obtained  from  these  as: 


Aco  =  — 

(38) 

td 

T 

N  = 

(39) 

2A  t 

The  velocity  and  flux  boundary  conditions  along  a  wave  generation  line  perpen¬ 
dicular  to  the  x-axis  may  be  obtained  from  the  surface  elevation  using  the  linear 
transfer  function  approach: 

N 

ua  (xg  >0  =  '5LTu  (<°,  )aj  cos  0  cos  [kj -(xg-x,)-  CO  jt  +  £j  ]  (40) 

j=i 
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Uf(xg,t)  =  [h  +  T\(xg,t)]ua(xg,t)  +  h3 


l 

6 


i^axx  (*,,0 +%,(*,,*)) 


(41) 


where 


fl  (*# » 0  =  2  «/  cos  [  Ay.  •  (xg  -  )  -  co/ +e.]  (42) 

7=1 

N 

uoxx  (Xg  ,t)  =  “X  a/y  COS3  0  COS  [*y  •  (xf  -  X,  )  -  (0/  +  8,  ]  (43) 

M 

N 

VwcyiXg’t)  =  Sin 2  0COS 0 C0S [Ay  ®/  +  Ey  ]  (44) 

7=1 

For  highly  nonlinear  irregular  waves  in  shallow  water1,  the  velocity  boundary 
conditions  would  have  to  be  modified  to  take  into  account  the  presence  of  lower 
and  higher  frequency  wave  components  induced  by  nonlinear  interactions 
between  the  primary  wave  components.  Second-order  Boussinesq  theory  can  be 
used  to  generate  the  velocities  and  flux  densities  associated  with  the  bound 
second-order  waves  along  the  wave  generation  boundary.  However,  this  theory  is 
valid  over  a  limited  range  of  wave  steepnesses  (H/L)  and  relative  depths  (h/L) 
because  it  only  includes  second-harmonic  terms.  Second-order  Stokes  theory,  for 
example,  cannot  accurately  describe  the  shape  of  cnoidal-type  waves  in  shallow 
water  when  the  Ursell  parameter  (HL2/h3)  is  large. 

A  different  approach  to  generating  steep  irregular  waves  in  shallow  water  is 
the  nonlinear  Fourier  method  of  Osborne  (1997),  in  which  an  irregular  wave  train 
is  represented  as  a  superposition  of  nonlinear  cnoidal  or  solitary  type  waves.  This 
approach,  however,  has  not  yet  been  implemented  in  BOUSS-2D. 

Irregular  Multidirectional  Waves.  Naturally  occurring  ocean  waves  exhibit 
a  pattern  that  varies  randomly  not  only  in  time  but  also  in  space.  The  wave 
energy  is  distributed  over  both  frequency  and  direction  and  can  be  described  in 
terms  of  a  directional  wave  spectrum,  S/to, 0),  which  is  the  product  of  the  fre¬ 
quency  spectrum,  S/co),  and  a  directional  spreading  function  Z>(co,0): 

S„(  ®,0)  =  ^(co)Z>(®,9)  (45) 


The  directional  spreading  function  is  non-negative  and  should  satisfy  the  follow¬ 
ing  relation: 


r 

J-n 


(46) 


1  A  useful  guideline  is  Hm JLp  >  0.025  tanh  kph,  where  is  the  significant  wave  height,  and  Lp 
and  kp  are  respectively  the  wavelength  and  wave  number  based  on  the  peak  frequency  of  the 
spectrum. 
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One  of  the  most  commonly  used  models  for  the  directional  spreading  function  is 
the  cosine-power  function  defined  as: 


D(Q)  =  T^-+1L  cos^e-e,)  for  1 0-0,1  <11/2  (47) 

Vtcn>+i/2)  p  p 

where  0^  is  the  principal  direction  of  wave  propagation  and  T  is  the  gamma  func¬ 
tion.  The  parameter  s  is  an  index  describing  the  degree  of  directional  spreading 
with  s  — >  °°  representing  a  unidirectional  wave  field.  Figure  4  shows  a  plot  of  the 
of  the  cosine  power  spreading  function  for  different  values  of  the  spreading 
index,  s. 


Figure  4.  Cosine-power  spreading  function  for  different  values  of  the  spreading 
index  s 

A  number  of  other  models  have  been  proposed  for  the  directional  spreading 
function  including  the  normal,  circular  normal,  and  wrapped-normal  distribution 
(Borgman  1969).  A  more  intuitive  and  universal  parameter  to  describe  the  degree 
of  directional  spreading  in  a  multidirectional  sea  state  is  the  standard  deviation  of 
the  directional  spreading  function,  oe,  which  is  defined  as: 


a 


2 

e 


J>0o+7T/2  ~ 


(48) 
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Given  the  same  standard  deviation,  there  is  very  little  difference  between  the 
shapes  of  different  proposed  forms  of  the  directional  spreading  functions  as 
shown  in  Appendix  C. 

Time-histories  of  the  velocity  boundary  conditions  can  be  synthesized  from 
the  directional  wave  spectrum  using  an  extension  of  the  linear  superposition 
technique  used  for  irregular,  unidirectional  waves.  The  directional  wave  spec¬ 
trum  is  initially  divided  into  N  frequency  bands  with  uniform  spacing,  Aco,  and  M 
direction  bands  with  uniform  spacing  A0.  The  water-surface  elevation  T|(jc,,f)  can 
be  represented  as  a  linear  superposition  of  periodic  wave  components  with  differ¬ 
ent  amplitudes  and  frequencies,  propagating  in  different  directions,  i.e., 

N  M 

T|(x, 0  =  SS av  cos | A  •  *  -  +  *y  ]  (49) 

i=\  j=l 

where  ky  =  (k,cos0y,  &;sin0y).  The  wave  amplitudes,  ay,  are  obtained  determin¬ 
istically  from  the  directional  wave  spectrum  as: 

Qy  =  ^2Sn  (cof ,  0y )  AcoA0  (50) 


while  the  phase  angles,  e,y,  are  randomly  selected  from  a  vmiform  distribution 
between  0  and  2k.  For  a  given  time  series  duration,  TD,  and  time-step,  At,  the 
values  of  N  and  Aco  can  be  obtained  in  the  same  manner  as  for  unidirectional 
waves  using  Equations  38  and  39  respectively. 

For  finite  duration  time  records,  the  double  summation  model  (Equation  49) 
produces  a  wave  field  that  is  spatially  nonhomogenous.  As  pointed  out  by 
Jeffreys  (1987),  this  is  because  the  phase  difference  between  wave  components 
with  identical  frequencies  but  propagating  in  different  directions  is  no  longer 
random  but  locked.  An  extreme  example  of  this  phenomenon  is  the  standing 
wave  pattern  that  is  produced  when  incident  waves  are  fully  reflected  by  a  verti¬ 
cal  wall.  The  wave  field  consists  of  two  waves  with  the  same  amplitude  and 
frequency  propagating  in  opposite  directions  with  a  constant  phase  difference 
between  them.  The  wave  heights  vary  from  zero  at  the  nodes  to  twice  the  inci¬ 
dent  wave  height  at  the  antinodes.  This  variability  is  obviously  reduced  by 
including  more  frequency  and  direction  components.  However,  as  pointed  by 
Stansberg  (1987)  and  others,  the  number  of  wave  components  required  to  reduce 
the  spatial  variability  in  wave  energy  to  within  acceptable  levels  is  quite  large 
[0(10,000)]. 

A  different  approach  to  producing  a  spatially  homogenous  multidirectional 
wave  field  is  to  assume  that  each  wave  component  has  a  unique  direction  of 
propagation.  Similar  to  an  irregular  unidirectional  wave  train,  the  water-surface 
elevation  is  given  by  Equation  36  with  kj  =  (kycos0„  £ysin0y).  Miles  (1989) 
discusses  four  different  methods  of  assigning  directions  to  each  frequency 
component.  In  BOUSS-2D,  the  random  direction  method  was  adapted  in  which 
wave  directions  are  selected  at  random  from  the  cumulative  distribution  of  the 
directional  spreading  function.  With  this  technique,  the  frequency  spectrum  is 
matched  exactly.  However,  the  directional  spreading  function  is  matched  only  in 
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a  global  sense  over  the  entire  record.  The  method  also  does  not  guarantee  a  sym¬ 
metric  directional  distribution  even  when  the  target  distribution  is  symmetric. 

The  velocity  and  flux  boundary  conditions  along  the  wave  generation  lines 
are  obtained  from  the  surface  elevation  time-histories  using  the  linear  transfer 
function  approach  described  previously  in  Equations  40  to  44. 


Internal  wave  generation  boundaries 

In  applications  where  there  is  significant  wave  reflection  from  bathymetric 
features  or  structures  within  the  computational  domain,  it  is  desirable  to  absorb 
the  waves  that  propagate  back  to  the  wave  generation  boundary  to  prevent  a 
buildup  of  wave  energy  inside  the  domain.  This  can  be  achieved  by  modifying 
the  boundary  conditions  along  the  wave  generation  boundary  to  simultaneously 
generate  and  absorb  reflected  waves  (e.g.,  Van  Dongeren  and  Svendsen  1997).  A 
different  approach  proposed  by  Larsen  and  Dancy  (1983)  is  to  generate  the 
waves  inside  the  computational  domain  and  absorb  reflected  waves  in  a  damping 
layer  placed  behind  the  generation  boundary.  This  approach  has  been  adopted  in 
BOUSS-2D  with  the  governing  equations  modified  to  allow  for  the  generation  of 
waves  inside  the  computational  domain. 

Consider  the  generation  of  waves  along  a  horizontal  line  by  a  distribution  of 
sources  that  extend  from  the  seabed  (z  =  - h )  to  the  free  surface  (z  =  T|).  The  veloc¬ 
ity  potential  associated  with  the  fluid  motion  satisfies  the  Laplace  equation 
everywhere  in  the  fluid  except  for  generation  line  (x  -  xg)  where  there  is  a  point 
source  of  fluid  mass.  The  governing  equation  for  the  fluid  motion  can  thus  be 
written  as: 


V2$  =  q(y,z,  t)  8(x -xg)  (51) 

where  q(y,z,t )  is  the  volume  flux  density.  Assuming  that  the  water  depth  is 
constant  along  the  generation  line,  a  modified  form  of  the  second-order  Taylor 
series  expansion  of  the  velocity  potential  about  an  arbitrary  elevation  z  =  za  in  the 
water  column  (Equation  1)  can  be  written  as: 

<t>  =  <t>«  +  ^[(za+h)2  ~(z  +  h)2]  [V2<|>a-g6(x-xg)]  (52) 

The  horizontal  fluid  velocities  are  obtained  from  the  velocity  potential  as: 

«  =  v<|>0  +  I[(za+/02-(z  +  /02]  V[V>a-*S <x-xg)]  (53) 

On  a  rectangular  grid  with  a  finite  grid  spacing  Ax,  the  delta  function  can  be 
replaced  with  2/Ax.  To  generate  waves  with  a  given  velocity  profile  u0(y,z,t ),  the 
mass  and  momentum  equations  along  the  grid  generation  line  and  adjacent 
velocity  points  can  be  written  as: 
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T|  +  V  •  P  u  dz  =  —  fn  u  dz 
'  J-*  A xJ-h 


(54) 


"c  +  £Vll  +  («a  V)  “a  +  ^[(Za+A)2-^]  [V  (V  ‘  "a,  “  2Uoa, , 1 **)]  =  0  (55) 


Some  previous  investigators  (e.g.,  Larsen  and  Dancy  1983)  introduced  a  point 
source  for  fluid  flow  into  the  continuity  equation.  It  is,  however,  important  to 
apply  a  correction  to  the  momentum  equations  to  account  for  higher-order  spatial 
derivatives  across  the  generation  line. 

Damping  regions 

Waves  propagating  out  of  the  computational  domain  are  absorbed  in  damp¬ 
ing  regions  placed  around  the  perimeter  of  the  computational  domain.  Damping 
layers  can  also  be  used  to  model  the  partial  reflection  from  harbor  structures 
inside  the  computational  area.  Artificial  dissipation  of  wave  energy  in  damping 
layers  is  achieved  through  the  introduction  of  a  term  proportional  to  the  surface 
elevation  into  the  right-hand  side  of  the  mass  equation: 

Fdamp,ng^  =  (56) 

and  a  term  proportional  horizontal  velocity  into  the  right-hand  side  of  the 
momentum  equation: 

F damping^  =  ~H(*)  »a  (57) 

where  |i(x)  is  the  damping  strength  with  units  of  s'1.  Introduction  of  the  damping 
terms  lead  to  the  nonconservation  of  mass  and  momentum  in  the  damping 
regions.  However,  extensive  tests  with  different  combinations  of  the  damping 
terms  showed  that  the  inclusion  of  damping  terms  in  both  the  mass  and  momen¬ 
tum  equation  was  more  effective  than  just  a  damping  term  in  the  momentum 
equation. 

Numerical  simulations  to  evaluate  the  performance  of  the  damping  layer 
showed  that  waves  could  be  effectively  damped  out  in  a  layer  half  a  wavelength 
wide  by  employing  a  quadratic  variation  of  p(jc)  with  a  peak  value  of  30 IT,  where 
T  is  the  wave  period.  The  damping  strength  has  thus  been  nondimensionalized  by 
30 IT,  i.e., 


M*)  =  ^(*)  (58) 

where  |X„/jc)  is  a  nondimensional  damping  coefficient  that  is  allowed  to  vary 
from  0  to  1.  Figure  5  shows  the  variation  of  the  reflection  coefficient  with  damp¬ 
ing  coefficient  for  different  relative  widths  (w/L)  of  the  damping  layer,  where  w 
is  the  width  of  the  damping  layer  and  L  is  the  wavelength.  It  should  be  pointed 
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Figure  5.  Variation  of  reflection  coefficient  with  damping  coefficient 

out  that  the  results  presented  in  Figure  5  are  for  normally  incident  waves.  Differ¬ 
ent  reflection  coefficients  would  be  obtained  for  obliquely  incident  waves. 


Flow  through  porous  structures 

The  Boussinesq  equations  have  also  been  modified  to  simulate  partial  wave 
reflection  and  transmission  through  surface-piercing  porous  structures  such  as 
breakwaters.  The  modified  equations  for  the  porous  region  can  be  written  in 
terms  of  either  the  velocity  in  the  pores  or  the  volume-averaged  (discharge) 
velocity.  In  order  to  easily  match  velocities  and  fluxes  across  the  water/porous 
structure  interface,  the  equations  for  the  porous  region  are  written  in  terms  of  the 
discharge  velocity.  Ignoring  inertial  effects,  the  modified  form  of  the  equations 
for  the  porous  region  are  obtained  by  replacing  u  with  u/n,  where  n  is  the  poros¬ 
ity,  and  including  a  term  to  account  for  energy  dissipation  inside  the  structure: 
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(59) 
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(60) 


««,,  +  «^VT1  +  «„V^j  +  + 

+  ^[{za+hf- fc2]v(V  uat)  +  nftua  +  nftua |»„ |  =  0 


where/  and /  are  laminar  and  turbulent  friction  factors  respectively.  Engelund 
(1953)  recommended  the  following  empirical  relationships  for  the  laminar  and 
turbulent  friction  factors: 
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(l-«)3  V 

K2  Z 


(61) 
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s  0-”) 1 

Po  n3  d 


(62) 


where  v  is  the  kinematic  viscosity  of  water,  d  is  the  characteristic  stone  size,  and 
aa  and  are  empirical  constants  that  range  from  780  to  1,500,  and  1.8  to  3.6 
respectively. 


Simulation  of  Wave  Runup 

The  runup  of  waves  on  shorelines  provides  an  important  boundary  condition 
for  predicting  wave-induced  currents  and  sediment  transport  in  the  surf  zone.  The 
runup  limit  is  also  important  for  determining  the  minimum  crest  elevation  of 
coastal  structures  to  prevent  overtopping  and/or  flooding.  A  simple  runup  scheme 
has  been  implemented  in  BOUSS-2D.  Dry  computational  cells  (land  points)  are 
assumed  to  be  porous  regions  where  the  phreatic  surface  elevation  and  volume- 
averaged  velocities  are  calculated  simultaneously  with  the  fluid  motion  in  the  wet 
cells.  When  the  phreatic  surface  elevation  exceeds  the  elevation  of  the  land  point 
by  a  specified  threshold,  the  porous  cell  is  considered  flooded  and  treated  as  a 
wet  cell  during  the  next  time-step.  Alternatively,  when  the  free  surface  elevation 
drops  below  a  specified  threshold  above  the  bottom  elevation  of  a  wet  cell,  the 
wet  cell  is  assumed  to  be  dry  and  treated  as  a  porous  cell  during  the  next 
time-step. 


Subgrid  Turbulence 

BOUSS-2D  optionally  provides  a  mechanism  to  simulate  the  turbulence  and 
mixing  that  occurs  in  regions  with  large  gradients  in  the  horizontal  velocities 
such  as  around  the  tips  of  breakwaters.  The  dissipation  term  is  identical  to  that 
used  for  wave  breaking,  i.e.,  Equation  16.  The  eddy  viscosity  is  given  by 
Smagorinsky’s  (1963)  formulation  with  the  turbulent  length  scale  proportional  to 
the  grid  size.  It  can  be  written  as: 
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V,  =  C'AxAy 


(j^a.x  )  (^a,.y  )  2  (^“  J ’  ^<*,x  ) 


where  Cs  is  an  empirical  constant. 


(63) 


26 


Chapter  3  Numerical  Solution 


4  Setting  Up  and  Running 
B0USS-2D 


B0USS-2D  is  a  high-resolution  wave  model  designed  for  investigating  com¬ 
plex  wave  transformation  problems  over  small  regions  (1-5  km).  It  is  ideally 
suited  for  applications  where  reflection,  diffraction,  and/or  nonlinear  interactions 
are  significant  such  as  near  coastal  inlets  and  harbors.  For  wave  propagation  over 
large  open  areas  where  the  processes  of  wind-wave  generation,  shoaling  and 
refraction  are  dominant,  phase-averaged  models  such  as  STWAVE  (Smith, 
Sherlock,  and  Resio  2001)  should  be  used.  For  periodic  harbor  oscillations, 
frequency-domain  models  such  as  the  elliptic  mild-slope  model  CGWAVE 
(Demirbilek  and  Panchang  1998)  should  be  used  since  time-domain  models  take 
longer  than  frequency-domain  models  to  attain  steady-state  conditions. 


Overview  of  Model  Setup 

The  basic  steps  involved  in  setting  up  and  running  BOUSS-2D  are: 

a.  Collect  ancillary  data  such  as  bathymetric  data  and  wave  climate 
information. 

b.  Select  portion  of  the  ocean  to  be  covered  by  the  numerical  model  and 
prepare  2-D  rectangular  grids  containing  information  on  the: 

( 1 )  Seabed  elevations  over  the  computational  area. 

(2)  Damping  values  for  the  absorption  of  outgoing  waves  at  model 
boundaries  (optional). 

(3)  Porosity  values  if  porous  structures  such  as  breakwaters  are  inside 
the  computational  domain  (optional). 

(4)  Tidal  current  distribution  over  the  computational  area  (optional). 

c.  Create  a  simulation  parameter  file  that  contains  all  the  information 
required  to  run  the  model.  This  can  be  done  using  the  DOS-based 
interactive  program  Pre-BOUSS2D  or  other  Windows-based  programs. 

d.  Run  the  BOUSS-2D  model. 

e.  Analyze  the  model  output  and  plot  the  results. 
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Collection  of  Bathymetric  and  Wave  Climate  Data 

The  first  step  in  the  modeling  process  is  the  collection  of  information  on  the 
seabed  elevations  over  the  computational  area.  This  can  be  obtained  from  hydro- 
graphic  surveys  of  the  area,  digitizing  nautical  charts,  or  from  digital  topographic 
databases  such  as  those  maintained  by  the  National  Oceanic  and  Atmospheric 
Administration  (NOAA). 

Wave  climate  information  close  to  the  project  site  can  be  obtained  from: 

a.  Long-term  wave  measurements  at  the  project  location. 

b.  Long-term  wave  measurements  in  deeper  water  with  buoys  such  as  those 
operated  by  NOAA. 

c.  Hindcast  from  long-term  wind  observations  using  models  such  as  Wave 
Model  (WAM),  Wave  Information  Steady  Wave  Model  (WISWAVE),  or 
Steady-State  Wave  Model  (STWAVE). 

Depending  on  the  source  of  the  wave  climate  data,  simpler  wave  propagation 
models  such  as  STWAVE  (Smith,  Sherlock,  and  Resio  2001)  or  Refraction/ 
Diffraction  Model  for  Spectral  Wave  Conditioning  (REFDIF-S)  (Kirby  and 
Ozkan  1994)  might  have  to  be  run  to  transfer  the  wave  climate  data  to  the 
boundary  of  the  computational  grid.  Statistical  techniques  (e.g.,  Borgman  1972) 
can  then  be  used  to  reduce  the  wave  climate  data  to  design  wave  conditions  with 
associated  return  periods. 

Information  on  tidal  water  levels  and  currents  near  the  project  site  can  be 
obtained  using  prediction  techniques  based  on  long-term  tidal  observations  (e.g., 
NOAA)  or  by  running  more  sophisticated  circulation  models  such  as  Advanced 
Circulation  (ADCIRC)  Model  for  Oceanic,  Coastal  and  Estuarine  Waters.  If  the 
current  speeds  exceed  10  percent  of  the  phase  velocity  of  the  waves,  the  currents 
could  significantly  modify  the  wave  field  and  a  2-D  circulation  model  would 
have  to  be  run  to  provide  the  spatial  distribution  of  current  speeds  over  the 
computational  area. 


Preparation  of  Bathymetric  Grid  File 

BOUSS-2D  simulates  wave  propagation  over  a  2-D  Cartesian  grid  as  shown 
in  Figure  6.  The  computational  domain  is  defined  as  a  rectangular  region  from 
(^origin,  y origin)  to  Cw,  >w)  with  uniform  grid  spacings  Ax  and  Ay  in  the  x  and  y 
directions,  respectively.  The  bathymetric  grid  represents  the  seabed  elevation  at 
each  node  of  the  grid  with  land  points  defined  as  positive  while  water  points  are 
defined  as  negative.  To  prepare  the  grid,  the  spatial  extent  and  spacing  of  the 
computational  grid  have  to  be  selected.  Factors  to  consider  in  the  selection  of  the 
grid  boundaries  and  spacing  include: 
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Figure  6.  Definition  sketch  for  computational  grid 

a.  The  grid  axes  should  be  aligned  as  much  as  possible  with  the 
predominant  wave  direction. 

b.  Points  of  interest  in  the  computational  domain  should  be  kept  at  least  one 
wavelength  away  from  the  boundaries  to  minimize  the  effect  of 
diffraction  into  the  damping  layers. 

c.  The  water  depth  should  be  uniform  along  the  wave  generation  boundary. 

d.  The  grid  spacing  should  be  chosen  to  resolve  the  shortest  wave  period  of 
interest  (Tmin)  in  the  shallowest  part  of  the  domain  (hmiD).  This  typically 
corresponds  to  having  at  least  eight  grid  points  per  wavelength  in  the 
shallow  regions  and  20  to  30  points  per  wavelength  at  the  peak  wave 
period  (Tp)  in  the  deep  regions  (/w)- 

e.  The  maximum  water  depth  has  to  be  less  than  half  the  shortest  wave¬ 
length  of  interest  (h^  <  L(Tmin)/2).  In  regions  with  relatively  deep  water 
along  the  offshore  boundary,  an  artificial  maximum  depth  can  be 
imposed  in  the  numerical  model  to  ensure  that  the  shortest  waves  of 
interest  are  resolved.  For  example,  if  we  are  interested  in  modeling  a  sea 
state  with  Tp  =10  s  and  =  6s  over  an  area  where  the  water  depth 
ranges  from  5  m  to  100  m,  a  maximum  depth  of  28  m  can  be  imposed 
corresponding  to  L(Tm -^12. 

f.  The  maximum  size  of  the  computational  area  should  be  based  on  the 
amount  of  available  computational  resources.  For  typical  PC’s  or  work¬ 
stations  with  100  to  200  MB  of  virtual  memory,  the  grid  size  should  not 
exceed  500  x  500  grid  points.  Larger  grids  can  be  used  at  the  cost  of 
much  longer  computational  times. 
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Once  the  grid  spacing  and  extent  have  been  chosen,  different  algorithms 
(linear,  inverse  distance,  kriging,  etc.)  can  be  used  to  interpolate  the  seabed 
elevation  data  (x,y,z)  onto  a  uniform  rectangular  grid.  Land  regions  and 
impermeable  harbor  structures  are  then  mapped  onto  the  grid  as  positive 
elevations.  The  grid  data  should  then  be  saved  as  an  ASCII  file  in  the  grid  file 
format  described  in  Appendix  D. 


Preparation  of  Damping  Grid  File 

The  external  boundaries  of  the  computational  grid  are  considered  to  be  fully 
reflecting  wall  boundaries  unless  it  is  specified  as  a  wave  generation  boundary. 
Damping  regions  have  to  be  placed  inside  the  domain  to  absorb  outgoing  waves. 
The  damping  regions  should  be  of  the  order  of  half  a  wavelength  wide  with  a 
quadratic  variation  of  the  nondimensional  damping  coefficient,  \in/x,y),  from  0 
to  1  at  the  wall  boundary. 

A  utility  program  GENDAMP  (Appendix  E)  has  been  provided  to  generate 
the  damping  file  given  the  location  of  the  damping  layers  (North,  South,  East,  or 
West  grid  boundaries),  the  width  of  the  damping  layer,  and  the  nondimensional 
damping  strength  at  the  end  wall.  The  damping  file  is  stored  as  an  ASCII  file  in 
the  grid  file  format  described  in  Appendix  D. 


Preparation  of  Porosity  Grid  File 

Porosity  layers  are  used  to  simulate  partial  reflection  and  transmission 
through  porous  structures  such  as  breakwaters.  The  porosity  grid  file  contains 
information  on  the  porosity  distribution  n(x,y)  over  the  computational  grid  with  a 
value  of  1 .0  for  water  points.  A  porosity  value  of  0.4  is  typically  used  for  break¬ 
waters.  Regions  with  small  porosity  values  {n  <  0.1)  such  as  breakwater  core 
layers  should  be  treated  as  impermeable  regions  and  mapped  as  land  points  in  the 
bathymetry  file. 

A  utility  program  MAP  POROSITY  (Appendix  E)  has  been  provided  to 
create  a  porosity  file  given  the  breakwater  boundaries  as  a  set  of  discrete  (x,y) 
points.  The  porosity  file  should  be  saved  as  an  ASCII  file  in  the  grid  file  format 
described  in  Appendix  D. 


Creation  of  Simulation  Parameter  File 

All  the  input  parameters  required  to  run  BOUSS-2D  is  stored  in  an  ASCII 
text  file,  referred  to  as  a  simulation  parameter  file  (.par).  The  file  is  a  free  format 
file  with  names  of  the  simulation  parameters  and  associated  values  not  con¬ 
strained  to  any  particular  rows  or  columns.  The  full  colon  (:)  character  is  used  at 
the  beginning  of  lines  containing  parameter  names  and  values,  and  the  pound  (#) 
character  is  used  to  delineate  sections  of  a  line  with  comments.  Due  to  the 
multitude  of  options  available  for  the  input  file,  it  is  recommended  that  the  inter¬ 
active  program  Pre_BOUSS2D  be  used  to  create  it.  Pre_BOUSS2D  can  be  run 
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from  an  MS-DOS  prompt  window  by  typing  the  program  name  including  the 
filepath,  e.g., 

C:\BOUSS2D\BIN\pre_bouss2d 

The  user  then  interactively  types  in  answers  on  the  keyboard  in  response  to 
questions  posed  by  the  program.  Default  answers  shown  in  square  brackets  [...] 
can  be  selected  by  hitting  the  enter  key. 

Step  1.  Simulation  Parameter  File 

•  Enter  the  name  of  the  output  simulation  parameter  file  (.par) 

Step  2.  Bathymetric  Information 

•  Enter  the  name  of  the  bathymetry  file  (.grd) 

•  Enter  storm  surge/tidal  water  level  offset  in  meters 

Step  3.  Damping  Information 

•  Enter  name  of  damping  file  (.grd)  if  there  are  damping  regions  in  the 
computational  domain 

Step  4.  Porous  Structure  Information 

Enter  the  following  information  if  there  is  a  porous  structure  in  the  computational 
domain 

•  Name  of  porosity  file  (.grd) 

•  Characteristic  stone  size  in  porous  layer  (m) 

Step  5.  Wave  Generation  Boundaries 

Waves  may  be  generated  along  one  or  two  boundaries.  The  boundaries  may  be 
external  (along  boundary  of  computational  grid)  or  internal  (within  the 
computational  grid).  Internal  wave  generation  boundaries  have  to  be  parallel  to 
one  of  the  grid  boundaries. 

For  each  wave  generation  boundary,  enter  the  following  information: 

•  Type  of  Wave  Generation  Boundary  (Extemal/Intemal) 

For  External  Wave  Generation  Boundary 

•  Select  grid  boundary  (North,  East,  South,  or  West) 

•  If  generation  line  exists  over  part  of  the  boundary,  enter  the  coordinates 
of  the  start  and  end  points. 

For  Internal  Wave  Generation  Boundary 

•  Enter  orientation  of  wave  generation  line  (East-West  or  North-South) 

•  Enter  the  x  or  y  coordinate  of  wave  generation  line 

•  If  generation  line  exists  over  part  of  boundary,  enter  the  coordinates  of 
the  start  and  end  points. 

Step  6.  Incident  Wave  Information 

Along  wave  generation  boundaries,  BOUSS-2D  requires  the  horizontal  velocities 
«o(0  and  flux  densities  «//)  normal  to  the  boundaiy  at  each  time-step.  The  time- 
histories  may  correspond  to  regular  or  irregular,  unidirectional  or  multidirectional 
waves.  The  velocity  and  flux  time  series  can  either  be  synthesized  from  para¬ 
metric  information  (wave  height,  period,  etc.)  or  derived  from  an  input 


Chapter  4  Setting  Up  and  Running  BOUSS-2D 


31 


surface-elevation  time  series.  The  information  required  to  generate  the  incident 
wave  boundary  conditions  is  described  as  follows: 

•  Wave  Synthesis  Options 

Synthesize  velocity  and  flux  time-histories  from  parametric  information 
Read  in  measured  surface-elevation  time  series.  The  input  file  has  to  be 
in  the  time  series  file  format  described  in  Appendix  D. 

•  Incident  Wave  Type 

Regular  Waves 

Irregular  Unidirectional  Waves 
Irregular  Multidirectional  Waves 

Regular  Waves 

•  Wave  Height  (m)  -  vertical  distance  between  the  wave  crest  and  trough. 
The  incident  wave  height  should  not  exceed  75  percent  of  the  breaking 
limit  given  by  the  Miche  criterion  as  Hbrcdk/L  =  0.14  tanh  kh.  When  the 
wave  height  exceeds  25  percent  of  the  breaking  limit,  nonlinear  effects 
become  important  and  the  time-histories  are  synthesized  using  the 
Boussinesq-Fourier  theory  discussed  in  Appendix  A. 

•  Wave  Period  (s)  -  time  interval  between  successive  wave  crests.  The 
incident  wave  period  has  to  be  greater  than  the  limit  imposed  by  the 
dispersive  properties  of  the  Boussinesq  equations,  i.e.,  the  wavelength 
calculated  from  the  linear  dispersion  relation  (Equation  10)  using  the 
wave  period  and  maximum  water  depth  has  to  be  greater  than  twice  the 
water  depth  (L  >  2/2max). 

•  Wave  Direction  (deg)  -  This  is  the  direction  the  waves  propagate  into 
the  computational  domain  from,  and  is  defined  in  a  clockwise  manner 
from  the  northern  boundary  of  the  grid  as  shown  in  Figure  6.  For 
example,  waves  propagating  into  the  domain  normal  to  the  western 
boundary  would  have  a  direction  of  270  deg.  The  incident  wave 
direction  also  has  to  be  within  ±  85  deg  of  the  normal  to  the  wave 
generation  boundary. 

•  Number  of  Wave  Cycles  -  This  determines  the  time  period  over  which 
output  parameters  such  as  the  significant  wave  height  and  mean  veloci¬ 
ties  are  calculated  after  steady-state  conditions  have  been  established  in 
the  computational  domain.  The  recommended  range  is  from  10  to  50. 

Irregular  Waves 

•  Type  of  Wave  Spectrum 

•  JONSWAP 

•  Bretschneider 

•  Pierson-Moskowitz 

•  TMA 

•  Ochi-Hubble 
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Different  spectral  parameters  will  be  input  depending  on  the  selected  wave 

spectrum.  The  most  commonly  used  parameters  are  described  in  the  follow¬ 
ing  paragraphs.  Detailed  information  on  the  different  wave  spectra  and  then- 

associated  parameters  are  provided  in  Appendix  B. 

•  Significant  Wave  Height,  Hmo  (m)  -  Defined  as  four  times  the  square 
root  of  the  zeroth  moment  of  the  spectrum,  i.e.,  Hmo=  4  yjm^ ,  where  mQ  = 

Jo  S^{f)df  .  Hmo  is  equivalent  to  the  average  height  of  the  highest  one- 

third  of  the  waves  (Hl/3)  for  moderate  sea  states  in  deep  water  with  a 
Rayleigh  wave  height  distribution.  The  input  significant  wave  height 
should  not  exceed  50  percent  of  the  breaking  limit  given  by  the  Miche 
criterion  based  on  the  peak  frequency  of  the  spectrum,  i.e.,  Hm0TmJLp  = 
0.07  tanh  kph. 

•  Peak  Wave  Period,  Tp  (s)  -  Inverse  of  the  cyclic  frequency  at  which  the 
wave  spectrum  is  a  maximum  ( Tp=\/fp ) 

•  Minimum  Wave  Period,  Tmin  (s)  -  The  incident  wave  spectrum  is  set  to 
zero  at  all  frequencies  greater  than  l/Tmin  with  the  truncated  wave  spec¬ 
trum  optionally  rescaled  to  match  the  target  significant  wave  height.  The 
minimum  wave  period  has  to  be  greater  than  the  limit  imposed  by  the 
dispersive  properties  of  the  Boussinesq  equations,  i.e.,  the  wavelength 
based  on  Tmin  and  has  to  be  greater  than  twice  the  maximum  water 
depth  [L  (7^)  >  2Amax]. 

•  Maximum  Wave  Period,  (s)  -  The  incident  wave  spectrum  is  set  to 
zero  at  all  frequencies  less  than  1/T^.  The  maximum  wave  period  is 
used  to  separate  infragravity  waves  from  wind-generated  waves  and  has 
a  default  value  of  25  s  (prototype). 

•  Wave  Direction  (deg)  -  This  is  the  direction  the  waves  propagate  into  the 
computational  domain  from,  and  is  defined  in  a  clockwise  manner  from 
the  northern  boundary  of  the  grid  as  shown  in  Figure  6.  For  example, 
waves  propagating  into  the  domain  normal  to  the  western  boundary 
would  have  a  direction  of  270  deg.  The  incident  wave  direction  also  has 
to  be  within  ±85  deg  of  the  normal  to  the  wave  generation  boundary. 

•  Synthesized  Time  Series  Duration  -  This  is  the  recycling  period  of  the 
Fourier  series  used  to  synthesize  the  time-histories  and  corresponds  to 
the  time  period  over  which  output  parameters  such  as  the  significant 
wave  height  and  mean  velocities  are  calculated  after  steady-state  condi¬ 
tions  have  been  established  in  the  computational  domain.  The  recom¬ 
mended  range  is  from  50  to  100  wave  periods. 

•  Random  Number  Seed  -  The  seed  of  the  random  number  generator  used 
to  select  phase  angles  for  the  different  Fourier  components.  Different 
seeds  can  be  used  to  obtain  different  time  records  from  the  same  wave 
spectrum. 
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Multidirectional  Waves 


In  addition  to  specifying  spectral  parameters,  the  user  has  to  select  additional 

parameters  to  describe  the  directional  distribution  of  wave  energy.  The 

different  directional  distributions  are  discussed  in  Appendix  C. 

•  Type  of  Directional  Spreading  Function 

Wrapped-Normal 

Cosine-Power 

•  Standard  Deviation,  ag  (deg)  -  The  standard  deviation  of  the  Wrapped- 
Normal  directional  spreading  function.  The  allowable  range  is  from 

5  deg  to  50  deg  with  10  deg  representing  narrow-band  swell-type  condi¬ 
tions  and  30  deg  representing  broader,  local  wind-generated  seas. 

•  Spreading  Index,  s  -  Describes  the  degree  of  directional  spreading  for  the 
Cosine-power  directional  distribution.  The  allowable  range  of  the 
spreading  index  is  from  0  to  65  with  s  =  2  representing  broad,  local 
wind-generated  seas  and  5  =  15  representing  narrow-band  swell-type 
conditions. 

•  Principal  Wave  Direction,  0P  (deg)  -  This  is  the  direction  corresponding 
to  peak  of  the  directional  spreading  function.  The  wave  direction  is 
defined  in  a  clockwise  manner  from  the  northern  boundary  of  the  grid  as 
shown  in  Figure  6.  The  principal  direction  has  to  be  within  ±  30  deg  of 
the  normal  to  the  wave  generation  boundary. 

•  Maximum  Propagation  Direction,  0™,*  (deg)  -  This  is  the  maximum 
wave  propagation  direction  relative  to  the  normal  to  the  wave  generation 
boundary  as  shown  in  Figure  7.  The  directional  distribution  is  truncated 
at  ±0max  and  renormalized.  As  pointed  out  by  Sand  and  Mynett  (1987), 
0max  also  defines  the  limited  area  in  the  computational  domain  over 
which  homogenous  conditions  exist,  i.e.,  all  wave  directions  are 
included.  Diffraction  effects  will  further  reduce  the  size  of  the  spatially 
homogenous  region.  The  recommended  value  is  two  to  three  times  the 
standard  deviation  of  the  directional  distribution. 

It  should  also  be  noted  that  when  using  more  than  one  wave  generation  boundary, 
different  spectral/directional  parameters  might  be  specified  for  the  individual 
boundaries.  However,  the  duration  of  the  synthesized  or  input  time  record  must 
be  the  same  for  both  boundaries. 

Step  7.  Current  Field 

If  tidal  currents  are  significant  over  the  computational  area,  the  spatial 
distribution  of  the  currents  U(x,y)  should  be  stored  as  an  ASCII  file  in  the  grid 
file  format  (.grd)  described  in  Appendix  D.  The  user  should  then  enter  the  name 
of  the  file  containing  the  velocity  data. 
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Figure  7.  Sketch  showing  spatially  homogenous  region  for  multidirectional 
waves 


Step  8.  Simulation  Parameters 

•  Duration  of  Numerical  Simulation  (s)  -  This  is  equal  to  time  required  for  the 
storm  waves  to  propagate  from  the  wave  generation  boundary  to  the  farthest 
points  of  interest  in  the  model  area  and  establish  steady-state  conditions,  i.e., 
warm-up  period  plus  the  duration  of  the  synthesized  or  input  time  record 
specified  in  Step  6. 

•  Simulation  Time-Step  (s)  -  The  time-step  should  be  selected  based  on  the 
stability  considerations,  i.e.,  the  Courant  number  Cr  has  to  be  less  than  1 .0: 


C  = 


c2  Ar 


l  l 

-+- 


Ax2  Ay2 


<  1.0 


where  C  is  the  phase  velocity  calculated  using  the  maximum  water  depth 
(hm ax)  and  the  average  zero-crossing  period  of  the  incident  waves.  It  is 
recommended  that  the  Courant  number  be  kept  within  the  range  0.5  to  0.7, 
which  typically  corresponds  to  30  to  50  points  per  wave  period.  For  plunging 
waves  on  steep  slopes,  it  may  be  necessary  to  use  up  to  100  points  per  wave 
period  to  simulate  the  rapid  changes  in  wave  shape  irrespective  of  the 
Courant  number. 


•  Chezy  Coefficient  -  Bottom  friction  factor  with  a  default  value  of  30.  Should 
be  kept  between  20  and  1,000. 

•  Model  Equations  Option  (Weakly  Nonlinear/Fully  Nonlinear)  -  Select  either 
weakly  nonlinear  model  of  Nwogu  (1993)  or  the  fully  nonlinear  model  of 
Nwogu  (1996).  The  fully  nonlinear  model  is  computationally  more  intensive 
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than  the  weakly  nonlinear  model  and  should  only  be  used  to  investigate 
highly  asymmetric  waves  in  shallow  water,  wave-induced  currents,  and 
wave-current  interaction. 

•  Wave  Breaking  Option  (Y/N)  -  Wave  breaking  should  be  enabled  if  the 
significant  wave  height  is  going  to  be  greater  than  half  the  water  depth  in  the 
shallowest  regions  of  the  computational  area  (Hm0  >  hmJ2). 

•  Turbulent  Length  Scale,  l,  (m)  -  Controls  the  rate  of  wave  energy  dissipation 
for  breaking  waves.  Should  be  set  equal  to  the  significant  wave  height  (/,  = 
Hm0)  for  irregular  waves  or  the  wave  height  {l,  -  H)  for  regular  waves. 

•  Smagorinsky  Constant  Cs -  Eddy  viscosity  coefficient  for  subgrid  turbulence. 
To  avoid  excessive  dissipation  of  the  waves,  Cs  should  be  kept  between  0  and 
0.5.  The  default  value  is  0.0. 

•  Wave  Runup  Option  (Y/N)  -  The  runup  scheme  used  in  the  numerical  model 
is  designed  to  simulate  subcritical  flow  conditions  on  mild  slopes.  It  cannot 
resolve  details  of  supercritical  flow  on  steep  slopes.  The  runup  option  should 
only  be  used  with  the  fully  nonlinear  model  equation  option. 

•  Minimum  Flooding  Depth  -  Parameter  used  control  the  stability  of  runup 
computations.  Its  value  should  be  of  the  order  of  one-hundredth  of  the 
incident  wave  height  (HI  100). 


Step  9.  Output  Data 

BOUSS-2D  calculates  the  time-dependent  evolution  of  the  water-surface  eleva¬ 
tion  and  horizontal  velocities  over  a  rectangular  grid.  For  most  simulations,  it 
would  require  an  excessive  amount  of  disk  space  to  store  the  surface  elevation 
and  velocity  data  at  every  grid  point  for  all  the  simulation  time-steps.  To  mini¬ 
mize  disk  storage  requirements,  the  program  outputs  the  time-averaged  values  of 
the  variables  over  the  entire  grid  or  time-histories  of  the  variables  at  specified 
grid  points  in  the  computational  domain.  All  output  files  start  with  a  prefix  (xxx) 
specified  by  the  user.  The  different  output  file  options  are  given  in  the  following 
paragraphs: 


2-D  Spatial  Output 

•  Mean  water  level  ff(x,y )  distribution  over  the  computational  grid 

(xxxmwl.grd)  -  The  free  surface  fluctuations  are  averaged  over  the  duration 
of  the  synthesized  or  input  time  record,  tw,  specified  in  Step  6,  i.e., 

,y,t) 

Nifiw 

where  N  =  tJAt  and  ts  is  the  duration  of  the  numerical  simulation  specified  in 
Step  8.  The  user  should  make  sure  that  the  duration  of  the  numerical  simula¬ 
tion  is  long  enough  to  establish  steady-state  conditions  in  the  numerical  wave 
basin. 
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•  Significant  wave  height  distribution  Hmo(x,y)  over  the  entire  computational 
grid  (xxxHs.grd)  -  The  significant  wave  height  is  calculated  as  four  times 
standard  deviation  of  the  water-surface  elevation  over  the  duration  of  the 
synthesized  or  input  time  record,  i.e.. 


Hmo(x,y)  =  4 


1 

N- 1 


V[{x,y,t)  -  r| 2(x,y) 


For  regular  waves,  the  wave  height  can  be  obtained  from  the  significant  wave 

height  as  Hmo/\l 2  .  It  should  noted  that  there  could  be  differences  between  the 
significant  wave  height  estimated  from  the  standard  deviation  of  the  time 
record  and  those  obtained  using  zero-crossing  analysis  especially  for  highly 
asymmetric  waves  in  shallow  water. 

•  The  time-averaged  or  mean  current  ua  (x,y)  over  the  entire  computational 
domain  (xxxuv.grd). 

1  U 

ua(,x,y)  =  —\ua{x,y,t) 

Nct 

The  time-averaging  procedure  is  also  carried  out  over  the  duration  of  the 
synthesized  or  input  time  record. 


Time  Series  Output 

•  Time  series  of  the  water-surface  elevation  7j(t)  and  two  components  of  the 
horizontal  velocity  ujt),  v„(t)  at  specified  (x,y)  locations  in  the  grid.  The 
time  series  files  are  named  xxx_eta.tsl,  xxx  u.tsl,  and  xxx  v.tsl. 

Animation  Output 

•  Time-dependent  output  of  the  water  surface  elevation  T]{x,y,f)  and  horizontal 
velocities  ua(x,y,t)  over  a  specified  area  and  time  period.  The  user  specifies 
the  coordinates  of  the  lower  left-hand  comer  (xuy\),  upper  right-hand  comer 
(*2,^2),  and  grid  skip  interval  for  the  output  grid  area,  and  the  start  time,  end 
time,  and  time-step  for  the  output  time  period.  The  files  are  named  xxx.eta 
and  xxx.uv  for  the  surface  elevation  and  velocities,  respectively.  The  anima¬ 
tion  files  could  be  large  and  should  be  used  judiciously.  The  size  of  the  sur¬ 
face  elevation  output  file  in  bytes  can  be  calculated  as  4  x  (number  of  grid 
points)  x  (number  of  time-steps).  The  velocity  file  is  twice  the  size  of  the 
surface  elevation  file. 

Step  10. 

After  entering  all  the  information,  Pre-BOUSS2D  creates  a  simulation  parameter 

file.  The  simulation  parameter  file  is  an  ASCII  file  that  can  be  edited  using 
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standard  text  editors  (WordPad,  etc.).  Changes  can  be  made  to  the  wave  simula¬ 
tion  parameters  (e.g.,  wave  height,  period,  direction,  time-step,  etc.)  and  the  file 
saved  under  a  new  name.  A  sample  output  of  the  simulation  parameter  file  is 
shown  as  follows: 


############################################################### 

#  BOUSS2D  Simulation  Parameter  File:  barbers_irr_tl2 . wsp 

#  Written  By:  John.  E.  Hacker 

#  Creation  Date:  Fri,  Sep  15,  2000 

#  Creation  Time:  04:16  PM 


############################################################### 

# 


#  Bathymetric  Grid  Parameters 

:  BATHY__FILE  bathy_barbers  .  grd 

: TIDAL_OFFSET  0  #  metres 

# 

#  Damping  Parameters 

: DAMPING_FILE  damp_barbers . grd 

# 

#  WaveMaker  #1  Parameters 
: START  WAVEMAKER 


WM_P0S_X1 

600 

# 

metres 

WM_P0S_Y1 

0 

# 

metres 

WM_POS_X2 

600 

# 

metres 

WM_POS_Y2 

2500 

# 

metres 

WAVE_TYPE  Irreg_ 

uni  # 

Regular, 

Irreg_uni,  or  Irreg_multi 

TIME_SERIES_OPTION 

synthesi: 

ze  # 

file  or  synthesize 

WAVE_DIRECTION 

270 

# 

degrees 

T I ME_S ER I E S_DURAT I ON 

1800 

# 

seconds 

RANDOM_NUMBER_SEED 

25136827 

SPECTRAL_TYPE 

JONSWAP 

#  JONSWAP,  BRET,  PM,  TMA  or  OCHI 

JONSWAP__OPTION 

Hs 

# 

Hs,  Sigma,  or  Gamma 

WAVE_HE I GHT_S IG 

3 

# 

metres 

WAVE__PERI  OD_PEAK 

12 

# 

seconds 

WAVE__PERIOD_MIN 

8 

# 

seconds 

WAVE_PERI OD_MAX 

24 

# 

seconds 

GAMMA 

3.3 

PHILLIPS 

0.0081 

RESCALE  SPECTRUM 

NO 

# 

YES  or  NO 

:  END  WAVEMAKER 
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#  Simulation  Parameters 

# 


: DURATION 
: TIME_STEP 
:  RAM  P_DURAT  I  ON 
:  CHEZY__COEFF 
: NONL I NE AR_0 PT I ON 
: WAVE_RUNUP_0 PT I ON 
: WAVE_BRE AKI NG_0 PT I ON 
: TURB_LENGTH_S  CALE 
# 

#  Output  Parameters 
: FILE_NAME_PREFIX 

: HS_FILE 
: MEAN_UV_FILE 

# 


2800 

# 

seconds 

0.2 

# 

seconds 

12 

# 

seconds 

30 

STRONG 

# 

WEAK  or  STRONG 

NO 

# 

YES  or  NO 

YES 

# 

YES  or  NO 

3 

# 

metres 

irr_h3tl2 

irr_h3t!2__Hs .  r2s  #  Hs  output  file 
irr_h3t!2_mean_UV. r2v  #  Mean  UV  output  file 


#  Surface  Elevation  (ETA)  File  Parameters 


: SAVE_ETA 
:  ETA_FILE 

:  START_TIME 

:  END_TIME 

:  SAVEJTIMEJBTEP 

:  SAVE_FULL_GRID 

;  GRID_STEP 
: END_SAVE_ETA 
# 

: SAVE_TIMESERIES 
:  TS_X 

:  TS_Y 

: END_TIMESERIES 

# 


irr_h3tl2 . wse 
0  #  seconds 

1800  #  seconds 

1  #  seconds 

YES  #  YES  or  NO 

1 


1220  #  metres 

890  #  metres 


Running  BOUSS-2D 

After  creation  of  the  simulation  input  file,  BOUSS-2D  can  be  run  from  an 
MS-DOS  command  prompt  window  by  typing: 

C: \BOUSS2D\BIN\bouss2d  input_f ile_name 
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After  carrying  out  computations  for  10  time-steps,  the  program  displays  an 
estimate  of  the  run  time.  The  run  time  could  vary,  however,  depending  on  other 
processes  running  on  the  computer  at  the  same  time. 


Time  Series  Data  Analysis 

BOUSS-2D  outputs  the  significant  wave  height  and  mean  current  distribution 
over  the  entire  computational  grid.  For  some  applications,  it  might  be  necessary 
to  analyze  the  time  series  data  at  specific  grid  points  to  obtain  frequency  or  direc¬ 
tional  wave  spectra,  wave  height  statistics,  reflection  coefficients,  etc.  Hughes 
(1993)  extensively  discusses  time  and  frequency  domain  analysis  techniques  for 
time  series  data  in  a  laboratory  setting.  Similar  techniques  can  be  applied  to  time 
series  output  from  the  numerical  model. 

Standard  statistical  analysis  techniques  can  be  applied  to  the  surface  eleva¬ 
tion  or  velocity  time  records  to  obtain  the  mean,  standard  deviation,  skewness, 
and  kurtosis  of  the  time  records.  The  skewness  of  the  surface  elevation,  r\,  and  its 
time  derivative,  i}„  can  be  used  to  quantify  the  degree  of  wave  asymmetry  or 
nonlinearity  in  shallow  water.  Zero-crossing  analysis  techniques  can  be  used  to 
define  the  heights  and  periods  of  individual  wave  cycles  within  an  irregular  wave 
record.  These  can  be  analyzed  to  obtain  parameters  such  as  the  root-mean-square 
wave  height,  H^,  average  height  of  the  highest  one-third  of  the  waves,  HVi,  or 
the  average  wave  period,  Tm. 

Standard  spectral  analysis  techniques  can  be  used  to  obtain  the  surface  eleva¬ 
tion  or  velocity  spectra.  These  can  be  analyzed  to  obtain  parameters  such  as  the 
significant  wave  height,  Hmo,  and  spectral  peak  period,  Tp.  Different  techniques 
have  been  proposed  for  the  estimation  of  directional  wave  spectra  using  either  an 
array  of  wave  gauges  (e.g.,  Nwogu  1989)  or  collocated  measurements  of  the 
surface  elevation  and  horizontal  velocities  (e.g.,  Nwogu  et  al.  1987).  Collocated 
time  series  of  the  surface  elevation  and  horizontal  velocities  can  also  be  used  to 
estimate  the  reflection  coefficients  using  the  technique  described  in  Hughes 
(1993). 
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Model  Validation 


Wave  Propagation  through  a  Breakwater  Gap 

We  initially  evaluated  the  ability  of  BOUSS-2D  to  simulate  the  propagation 
of  waves  through  a  gap  into  a  rectangular  harbor  basin.  The  basin  is  1,200  m 
wide,  600  m  long,  and  10  m  deep.  The  width  of  the  opening  at  the  harbor 
entrance  is  120  m. 

The  first  case  considered  is  a  regular  wave  with  period  T=  7  s  propagating  in 
a  direction  normal  to  the  breakwaters.  The  corresponding  ratio  of  the  gap  width, 
B,  to  wavelength,  L,  is  2.  Boussinesq  model  simulations  were  carried  out  using 
grid  spacings  Ax  =  Ay  =  3  m  and  time-step  size  At  =  0. 1 5  s.  A  60-m-wide  damp¬ 
ing  layer  was  placed  around  the  perimeter  of  the  basin  to  absorb  outgoing  waves. 
Figure  8  shows  a  snapshot  of  the  instantaneous  water-surface  elevation  produced 
by  the  BOUSS-2D  model. 

The  normalized  wave  height  distribution  predicted  by  BOUSS-2D  is  com¬ 
pared  with  the  numerical  solution  of  Isaacson  and  Qu  (1990)  in  Figure  9. 

Isaacson  and  Qu  (1990)  used  a  boundary  integral  technique  to  solve  the  2-D 
Helmholtz  equation,  which  is  a  reduced  form  of  the  mild-slope  equation  for 
water  of  constant  depth.  Good  agreement  is  generally  observed  between  the  wave 
height  predictions  from  the  numerical  models.  Small  oscillations  can  be  seen  in 
the  Boussinesq  model  predictions,  especially  for  the  smaller  wave  height  con¬ 
tours.  This  is  due  to  partial  wave  reflection  from  the  radiation  boundaries.  Reflec¬ 
tion  coefficients  of  the  order  of  5  to  10  percent  at  the  boundaries  have  been 
observed  to  cause  such  oscillations  in  wave  height  contour  lines  (Isaacson  and 
Qu  1990) 

We  next  investigated  the  propagation  of  irregular  multidirectional  waves 
through  the  gap.  The  sea  state  is  characterized  by  a  JONSWAP  wave  spectrum 
with  peak  period  Tp  =  7  s  and  peak  enhancement  factor  y=  3.3.  A  wrapped- 
normal  distribution  with  a  standard  deviation  of  20  deg  was  used  for  the  direc¬ 
tional  distribution  of  wave  energy.  The  double-summation  method  was  used  to 
synthesize  time-histories  of  velocity  fluxes  along  the  incident  wave  boundary  for 
the  BOUSS-2D  simulations. 

Figure  10  shows  a  snapshot  of  the  instantaneous  water-surface  elevation 
produced  by  the  BOUSS-2D  model.  The  normalized  wave  height  distribution 
predicted  by  BOUSS-2D  is  compared  with  the  numerical  solution  of  Isaacson 
and  Qu  (1990)  in  Figure  1 1 .  Good  agreement  is  observed.  As  expected,  the 
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y/L  (m) 


Figure  8.  3-D  view  of  instantaneous  water-surface  elevation  for  regular  waves 

propagating  through  a  breakwater  gap  (7=  7  s,  h  =  10  m,  B/L  =  2) 


Figure  9.  Relative  wave  height  contours  for  regular  waves  propagating  through 
a  breakwater  gap  (7=  7  s,  h  =  10  m,  B/L  =  2) 
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Figure  10.  3-D  view  of  instantaneous  water-surface  elevation  for  multidirectional 
waves  propagating  through  a  breakwater  gap  (Tp  =  7  s,  oe  =  20°,  h  = 
10  m,  B/Lp-2) 


Figure  1 1 .  Relative  wave  height  contours  for  multidirectional  waves  propagating 
through  a  breakwater  gap  (Tp  =  7  s,  o6  =  20°,  h  =  10  m,  B/Lp  =  2) 
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directional  spreading  of  wave  energy  leads  to  larger  wave  heights  in  the  sheltered 
area  behind  the  breakwaters  and  smaller  wave  heights  along  the  principal  direc¬ 
tion  of  the  waves. 


Multidirectional  Wave  Propagation  over  a  Shoal 

Laboratory  experiments  on  the  transformation  of  irregular  multidirectional 
waves  over  an  elliptical  shoal  were  carried  out  by  Vincent  and  Briggs  (1989). 
The  experimental  layout  is  shown  in  Figure  12,  and  consists  of  an  elliptical  shoal 
placed  in  a  0.4572-m-deep  basin.  The  boundary  of  the  shoal  is  an  ellipse  defined 
by: 


3.96 


V 

/ 


+ 


IZlc 

3.05 


=  1 


(64) 


where  (xc,  yc )  are  the  coordinates  of  the  center  of  the  shoal  and  are  given  by  xc  = 
13.72  m  and yc  =  6.10  m.  The  water  depth  over  the  shoal  is  given  by: 


h(x,y )  =  0.9144  -  0.7620 


x-x. 


\2 


V  4-95  7 


y-yc 

3.81 


V 


,  0.5 


(65) 


The  minimum  water  depth  over  the  shoal  is  0.1524  m.  Tests  were  carried  out 
for  various  regular  and  irregular,  unidirectional  and  multidirectional  waves.  For 
irregular  waves,  the  TMA  spectrum  was  used  to  describe  the  frequency  distribu¬ 
tion  of  wave  energy,  while  the  wrapped-normal  distribution  was  used  for  the 
directional  spreading  function.  Water-surface  elevation  data  were  collected  at 
five  transects  located  at  distances  of  3.05  m,  6.10  m,  9.14  m,  12.19  m,  and 
15.24  m  from  the  wavemaker  as  shown  in  Figure  12. 

The  numerical  basin  for  the  BOUSS-2D  model  simulations  is  31.5  m  wide, 
27  m  long,  with  a  uniform  grid  spacing  of  0.1  m.  Two-meter-wide  damping 
layers  were  placed  around  the  perimeter  of  the  basin  to  absorb  outgoing  waves. 
Two  representative  test  cases  were  selected  for  the  model-data  comparisons.  Test 
case  N1  is  characterized  by  a  TMA  spectrum  with  significant  wave  height  Hmo  = 
0.0775  m,  peak  period  7^,  =  1.3  s,  peak  enhancement  factor  y  =  2,  and  a  narrow 
directional  distribution  with  standard  deviation  G(,  =  10  deg.  Test  case  B1  is 
characterized  by  a  TMA  spectrum  with  significant  wave  height  Hmo  =  0.0775  m, 
peak  period  Tp  =  1.3  s,  peak  enhancement  factor  y=  2,  and  a  broad  directional 
distribution  with  standard  deviation  a9  =  30  deg.  Time-histories  of  the  velocity 
boundary  conditions  along  the  wavemaker  were  synthesized  for  a  duration  of 
130  s,  corresponding  to  100  wave  periods.  Simulations  were  carried  out  using  a 
time-step  size  At  =  0.025  s. 
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Figure  12.  Plan  view  of  bathymetry  and  layout  for  Vincent-Briggs  shoal 
experiments 

Figure  1 3  shows  a  snapshot  of  the  instantaneous  water-surface  elevation  over 
the  shoal  predicted  by  the  BOUSS-2D  model  for  test  case  Nl.  The  focusing  of 
wave  energy  behind  the  shoal  can  clearly  be  seen.  The  corresponding  2-D  map  of 
the  normalized  wave  height  distribution  over  the  computational  grid  is  shown  in 
Figure  14.  The  predicted  wave  height  variation  along  transects  3  and  4  are  com¬ 
pared  to  the  experimental  results  in  Figures  15  and  16,  respectively.  BOUSS-2D 
significantly  overpredicted  the  wave  height  amplification  along  transect  3  with 
differences  of  the  order  of  20  to  30  percent.  A  better  match  was  obtained  along 
transect  4  where  the  wave  height  amplification  was  lower.  The  reason  for  such 
large  discrepancies  along  transect  3  is  still  unclear.  Similar  discrepancies  were 
obtained  with  the  CGWAVE  model  (Demirbilek  and  Panchang  1998)  as  shown 
in  Figures  15  and  16.  For  the  broader  directional  distribution  test  case  Bl,  rea¬ 
sonable  agreement  was  obtained  between  the  model  and  the  lab  data  as  shown  in 
Figures  17  and  18  with  differences  of  the  order  of  10  percent. 


Wave  Breaking  in  Bimodal  Sea  States 

Bimodal  sea  states  consisting  of  swell  and  local  wind-generated  components 
occur  frequently  along  most  U.S.  coastlines  (Thompson  1980).  Smith  and 
Vincent  (1992)  investigated  the  shoaling  and  breaking  of  dual-wave  systems  on  a 
constant  slope  beach.  They  found  that  the  higher  frequency  component  decayed 
much  faster  in  the  presence  of  the  low-frequency  waves,  while  the  lower  fre¬ 
quency  component  appeared  to  be  unaffected  by  the  presence  of  high-frequency 
waves.  Smith  and  Vincent  (1992)  investigated  the  applicability  of  spectral  energy 
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Figure  13.  3-D  view  of  multidirectional  wave  propagation  over  a  shoal  for  test 
case  N1  (Hmo  =  0.0775  m,  Tp=  1.3  s,  <re  =  10  deg) 


Figure  14.  Normalized  wave  height  distribution  for  multidirectional  wave 
propagation  over  a  shoal  for  test  case  N1 
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Figure  17.  Normalized  wave  height  distribution  along  transect  3  for  test  case  B1 


Figure  18.  Normalized  wave  height  distribution  along  transect  4  for  test  case  B1 
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conservation  models  with  a  depth-limited  breaking  criterion  and  simpler  tech¬ 
niques  such  as  shoaling  the  component  wave  systems  independently  of  each 
other  and  superimposing  the  results.  None  of  these  techniques  could  reproduce 
the  observed  changes  to  the  wave  spectrum. 

Boussinesq  model  simulations  were  carried  out  for  test  case  7  in  the  Smith 
and  Vincent  experiments.  The  sea  state  is  composed  of  two  irregular  wave 
components  described  by  TMA  spectra  with  Hm0il  =  7.6  cm,  Tp<\  =  2.5  s,  Hm0t 2  = 
13.2  cm,  TPii  =  1.75  s,  and  y  =  20.  The  numerical  wave  flume  is  41  m  long  and 
consists  of  a  20-m-long  constant-depth  region  and  a  1 :30  planar  beach.  The  water 
depth  in  the  constant  depth  section  of  the  flume  was  0.61  m.  The  measured  time 
series  at  the  offshore  gauge  was  used  to  generate  velocity  boundary  conditions 
for  the  numerical  model.  A  grid  spacing  of  0.125  m  and  time-step  size  of  0.04  s 
was  used  for  the  numerical  simulations. 

Figures  19  to  22  show  comparisons  of  the  measured  and  predicted  wave 
spectra  at  gauges  located  in  water  depths  of  61  cm,  18.3  cm,  9.1  cm,  and  6.1  cm. 
Although  the  higher  frequency  component  dominates  the  incident  wave  spectrum 
with  67  percent  of  the  wave  energy  at  h  =  61  cm,  its  energy  is  preferentially  dis¬ 
sipated  with  the  lower  frequency  component  becoming  dominant  at  the  shallow 
depth  of  6.1  cm.  The  Boussinesq  model  is  able  to  reproduce  the  observed  trends 
in  the  experimental  data.  The  primary  reason  for  the  preferential  reduction  in 
energy  of  the  high-frequency  component  is  due  to  the  nonlinear  cross-spectral 
transfer  of  energy  during  the  shoaling  process. 


Figure  1 9.  Measured  and  predicted  wave  spectra  at  Gauge  1  for  bimodal  sea 
state  shoaling  on  a  constant  slope  beach 
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Figure  20.  Measured  and  predicted  wave  spectra  at  Gauge  4  for  bimodal  sea 
state  shoaling  on  a  constant  slope  beach 


Figure  21.  Measured  and  predicted  wave  spectra  at  Gauge  7  for  bimodal  sea 
state  shoaling  on  a  constant  slope  beach 
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Figure  22.  Measured  and  predicted  wave  spectra  at  Gauge  9  for  bimodal  sea 
state  shoaling  on  a  constant  slope  beach 


Rip  Currents  on  Barred  Beaches 

Rip  currents  are  narrow,  jet-like  currents  that  flow  out  from  the  surf  zone 
towards  the  open  ocean.  Rip  currents  are  typically  generated  when  there  are 
alongshore  variations  in  the  wave  breaking  location.  Alongshore  variations  could 
be  associated  with  bathymetric  effects  (e.g.,  gaps  in  sandbars),  the  presence  of 
structures  (piers,  jetties,  etc.)  or  edge  waves,  generated  by  the  trapping  of 
reflected  waves  near  the  shoreline.  One  characteristic  feature  of  rip  currents  is 
their  unsteady  nature.  The  currents  tend  to  occur  episodically  as  well  as  oscillate 
in  the  horizontal  plane  (Smith  and  Largier  1995). 

Haller,  Dalrymple,  and  Svendsen  (1997)  carried  out  laboratory  experiments 
to  investigate  the  generation  of  rip  currents  on  a  barred  beach  with  rip  channels. 
The  bathymetry  consists  of  a  1 :30  constant  slope  beach  with  a  superimposed 
longshore  bar.  The  bar  has  two  1.8-m-wide  gaps  as  shown  in  Figure  23.  The 
water  depth  is  40  cm  at  the  seaward  boundary  and  5  cm  at  the  top  of  the  bar.  The 
Boussinesq  model  was  applied  to  a  regular  wave  test  case  with  height 
H=  4.5  cm,  and  period  T=  1  s.  Numerical  simulations  were  carried  out  using 
grid  spacings  Ax  =  Ay  =  0.06  m  and  time-step  size  At  =  0.015  s. 

Figure  24  shows  a  snapshot  of  the  instantaneous  water-surface  elevation  70  s 
into  the  simulation.  Waves  propagating  over  the  bar  break  on  the  bar  while  those 
propagating  through  the  gap  break  closer  to  the  shoreline.  This  sets  up  a  spatial 
variation  of  the  mean  water  level  and  drives  a  time-varying  circulation  pattern. 
Figure  25  shows  a  2-D  map  of  the  mean  currents  averaged  over  50  wave  cycles 
(from  t  =  150  s  to  t  =  200  s).  As  observed  in  the  laboratory  experiments,  there  are 
two  pairs  of  counterrotating  circulation  cells.  The  primary  circulation  consists  of 
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Figure  23.  Plan  view  of  bathymetry  for  rip  current  experiments 


Figure  24.  3-D  view  of  wave  propagation  over  a  barred  beach  with  a  rip  channel 
(H  =  5  cm,  7=  1  s) 
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Figure  25.  Time-averaged  rip  current  pattern  at  f  =  200  s 

opposing  longshore  currents  that  converge  and  flow  seawards  through  the  gap  in 
the  sandbar.  The  secondary  circulation  cell  occurs  closer  to  the  shoreline. 


Wave-Current  Interaction 

The  presence  of  currents  can  significantly  modify  the  wave  field  inside 
entrance  channels  and  harbors.  This  is  particularly  true  for  adverse  currents 
where  the  interaction  shortens  and  steepens  the  waves,  leading  to  potentially 
hazardous  navigation  conditions.  Smith  et  al.  (1998)  carried  out  laboratory 
experiments  to  investigate  wave-current  interaction  in  an  idealized  inlet.  The 
bathymetry  and  layout  used  for  the  numerical  simulations  is  shown  in  Figure  26. 
The  wave  basin  is  24  m  wide  and  18.25  m  long.  The  water  depth  variation  adja¬ 
cent  to  the  inlet  is  given  by  the  equilibrium  beach  profile: 

h  =  0.036 y067  +  0.022  (66) 

The  equilibrium  profile  extends  to  a  water  depth  of  0.205  m  beyond  which  it  is 
linearly  transitioned  at  a  slope  of  1 :7  to  the  constant  depth  region  of  0.325  m.  The 
entrance  channel  is  lined  by  two  parallel  jetties  that  extend  a  distance  of  3.8  m 
from  the  shoreline  (y  =  0)  with  an  entrance  width  of  4  m. 
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Figure  26.  Bathymetry  of  idealized  inlet  for  wave-current  interaction  study 


A  2-D  hydrodynamic  model  was  run  to  provide  the  current  field  for  use  in 
the  wave  model  simulations.  Figure  27  shows  a  map  of  the  current  field 
generated  by  the  model.  Boussinesq  model  simulations  were  carried  out  using 
grid  spacings  Ac  =  Ay  =  0.075  m  and  time-step  size  At  =  0.015  s.  The  incident 
wave  conditions  were  characterized  by  a  TMA  spectrum  with  Hm0  =  0.055  m, 
Tp=  1.4  s,  and  y=  3.3.  The  predicted  wave  height  distributions  for  current  speeds 
U=  0  m/s  and  U=  0.24  m/s  are  shown  in  Figures  28  and  29,  respectively.  The 
Boussinesq  model  predicted  an  increase  of  15  to  20  percent  in  the  significant 
wave  height  near  the  entrance  channel.  The  laboratory  measurements,  however, 
showed  a  slight  decrease  in  wave  height  through  the  channel  due  to  the  effect  of 
wave  breaking.  A  more  detailed  investigation  of  the  wave-breaking  criterion  and 
energy  dissipation  rates  in  the  presence  of  currents  will  be  carried  out  to  improve 
simulations  of  wave-current  interaction  in  the  Boussinesq  model. 


Wave  Transformation  Near  Ponce  de  Leon  Inlet, 
Florida 

Ponce  de  Leon  Inlet,  an  inlet  leading  into  the  Halifax  and  Indian  Rivers  in 
Florida  has  a  complex  bathymetry  featuring  a  large  ebb  shoal,  a  navigation 
channel,  and  a  jetty.  A  1:100  laboratory  model  study  of  a  4.2-km  by  1.4-km 
region  near  the  inlet  was  carried  out  at  the  U.S.  Army  Engineer  Research  and 
Development  Center  (ERDC)  Coastal  and  Hydraulics  Laboratory  (CHL).  Tests 
were  conducted  for  sea  states  with  different  wave  heights,  peak  periods,  spectral 
widths,  and  directional  spread.  The  TMA  spectrum  was  used  to  describe  the 
frequency  distribution  of  wave  energy,  while  the  wrapped-normal  distribution 
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Figure  28.  Predicted  wave  height  distribution  near  inlet  for  test  case  without 
currents  (Hm0  =  0.055  m,  Tp  =  1 .4  s,  U  =  0  m/s) 
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Figure  29.  Predicted  wave  height  distribution  near  inlet  for  test  case  with 
currents  (F/mo  =  0.055m,  Tp  =  1 .4  s,  U  =  0.24  m/s) 


was  used  for  the  directional  spreading  function.  Water-surface  elevation  data 
were  recorded  at  30  locations. 

BOUSS-2D  was  run  for  one  of  the  multidirectional  test  cases  characterized 
by  Hmo  =  0.95  m,  Tp  =  10  s,  y=  5,  and  directional  spread  <r0  =  20  deg  (Test 
No.  11).  The  bathymetry  used  for  the  numerical  simulations  is  shown  in  Fig¬ 
ure  30.  Boussinesq  model  simulations  were  carried  out  using  Ax  =  Ay  =  5  m  and 
At  =  0.15  s.  Time-histories  of  the  velocity  boundary  conditions  along  the  wave- 
maker  were  synthesized  using  the  double-summation  method  for  a  duration  of 
500  s,  corresponding  to  50  wave  periods. 

A  snapshot  of  the  instantaneous  water-surface  elevation  produced  by  the 
BOUSS-2D  model  is  shown  in  Figure  31.  The  corresponding  significant  wave 
height  distribution  is  shown  in  Figure  32.  Two  distinct  wave-focusing  regions 
can  be  observed  on  the  shoal.  The  predicted  wave  height  variations  along  the 
offshore  and  nearshore  gauge  arrays  are  compared  to  the  experimental  results  in 
Figures  33  and  34,  respectively.  BOUSS-2D  reasonably  reproduced  the  wave 
height  variation  along  the  offshore  array,  although  differences  of  the  order  of 
20  percent  exist  at  a  couple  of  gauge  locations.  For  the  nearshore  array,  excellent 
agreement  is  obtained  between  BOUSS-2D  model  predictions  and  the  experi¬ 
mental  data. 
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Figure  30.  Ponce  de  Leon  Inlet  model  bathymetry 


Figure  31 .  3-D  view  of  multidirectional  wave  propagation  near  Ponce  de  Leon 
Inlet  (Hmo  =  0.95  m,  Tp  =  10  s,  cre  =  20  deg) 
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Figure  32.  2-D  map  of  wave  height  distribution  predicted  by  Boussinesq  model 
(Hmo  =  0.95  m,  Tp  =  10  s,  cre  =  20  deg) 


Figure  33.  Measured  and  predicted  wave  height  distribution  along  the  offshore 
gauge  array  (Hmo  =  0.95  m,  Tp  =  10  s,  Oe  =  20  deg) 
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Figure  34.  Measured  and  predicted  wave  height  distribution  along  the  nearshore 
gauge  array  ( Hmo  =  0.95  m,  Tp  =  10  s,  <re  =  20  deg) 


Wave  Disturbance  in  Barbers  Point  Harbor,  Hawaii 

Barbers  Point  Harbor  is  a  small  harbor  in  Hawaii  that  has  experienced 
occasional  long-period  harbor  oscillation  problems.  The  main  harbor  basin  is 
1 1.6  m  deep,  and  is  connected  to  the  Pacific  Ocean  through  a  12.8-m-deep, 
1.3-km-long  entrance  channel.  The  water  depths  are  shallow  (2  to  10  m)  in  the 
vicinity  of  the  entrance  channel,  leading  to  the  nonlinear  generation  of  free  and 
forced  long-period  waves  near  the  harbor  entrance.  Free  long  waves  could  be 
resonantly  amplified  inside  the  harbor  basin  if  the  long-wave  periods  are  close  to 
any  of  the  natural  modes  of  oscillation  of  the  basin. 

BOUSS-2D  was  used  to  investigate  to  the  resonant  harbor  oscillation  periods 
and  the  wave  height  amplification  factors.  A  three-dimensional  view  of  the 
bathymetry  and  harbor  layout  used  for  the  simulations  is  shown  in  Figure  35.  The 
offshore  boundary  was  truncated  at  a  water  depth  of  50  m.  Numerical  simulations 
were  carried  for  regular  waves  with  periods  ranging  from  50  s  to  200  s.  Wave 
height  information  was  output  at  two  gauges  outside  the  harbor  and  four  gauges 
inside  the  harbor  as  shown  in  Figure  36. 
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The  predicted  wave  heights  at  two  gauge  locations  inside  the  basin,  normal¬ 
ized  by  the  incident  wave  height,  are  compared  with  predictions  from  the  elliptic 
mild-slope  model  CGWAVE  (Demirbilek  and  Panchang  1998)  in  Figures  37  and 
38.  The  harbor  exhibits  a  number  of  distinct  periods  of  oscillation  (60  s,  83  s, 

1 10  s,  130  s,  143  s,  and  190  s).  The  Helmholtz  or  pumping  mode  of  the  basin 
occurs  at  a  period  of  around  900  s.  Good  agreement  is  generally  observed 
between  wave  height  amplification  factors  predicted  by  both  models.  It  should  be 
pointed  out  that  it  takes  longer  to  attain  steady-state  conditions  in  time-domain 
models,  especially  for  resonant  oscillations.  Figure  39  shows  a  plot  of  the  time 
history  at  Gauge  5  for  one  of  the  resonance  periods  (T=  60  s).  Steady-state 
conditions  are  attained  approximately  30  wave  periods  after  the  waves  initially 
arrived  at  the  gauge  location. 

Although  linear,  frequency-domain  models  are  computationally  more  effi¬ 
cient  at  predicting  harbor  resonance  periods  and  amplification  factors,  they 
cannot  predict  the  magnitude  of  the  long-period  wave  energy  inside  a  harbor 
from  a  given  offshore  wind-wave  spectrum.  To  overcome  this  deficiency, 
Okihiro,  Guza,  and  Seymour  (1993)  used  an  ad-hoc  coupling  of  a  nonlinear 
model  for  the  generation  of  bound  long  waves  outside  a  harbor  with  a  linear 
model  for  the  amplification  of  long  waves  inside  the  harbor.  The  complex  nature 
of  bathymetry  outside  Barbers  Point  Harbor  makes  it  difficult  to  quantify  the 
relative  amount  of  long-wave  energy  outside  the  harbor  that  is  freely  propagating 
into  the  harbor.  The  entrance  channel  is  much  deeper  than  the  surrounding  areas. 
Free  long  waves  would  be  generated  along  the  steep  side  slopes  of  the  entrance 
channel  as  well  as  reflected  from  shoreline.  The  long  waves  would  thus  be  propa¬ 
gating  over  a  wide  range  of  directions. 

We  investigated  the  ability  of  the  Boussinesq  model  to  simultaneously  model 
the  nonlinear  generation  of  long  waves  by  storm  waves  propagating  from  deep  to 
shallow  water,  the  diffraction  of  both  short  and  long  period  waves  into  the  harbor 
basin,  and  the  resonant  amplification  of  long  waves  inside  a  harbor.  We  initially 
considered  a  bichromatic  wave  train  with  component  periods  7\  =  12  s,  T2  = 

1 3 .46  s,  and  heights  H\=H1=L  1.5  m.  The  group  period  of  1 1 0  s  corresponds  to 
one  of  the  natural  periods  of  oscillation  of  the  basin.  Numerical  simulations  were 
carried  out  with  Ax  =  Ay  =  10  m  and  A t  =  0.2  s.  The  simulated  surface  elevation 
time-histories  at  the  offshore  Gauge  1,  harbor  entrance  Gauge  2,  and  harbor  basin 
Gauges  3  and  5  are  shown  in  Figure  40.  The  long-period  component,  obtained  by 
applying  a  low-pass  filter  (T>  25  s),  is  also  shown  in  the  figures.  It  can  be  seen 
that  nonlinear  interactions  during  the  shoaling  process  lead  to  an  amplification  of 
the  long-period  wave  component  between  the  offshore  Gauge  1  which  is  in  50  m 
of  water,  and  the  harbor  entrance  Gauge  2  which  is  7  m  of  water.  Inside  the 
harbor  basin,  the  long  waves  are  further  amplified  and  dominate  the  harbor 
response  at  Gauge  5. 

Although  bichromatic  waves  are  useful  for  demonstrating  the  importance  of 
nonlinear  wave-wave  interactions  in  harbor  response,  natural  sea  states  are 
irregular  with  wave  energy  distributed  over  a  large  number  of  frequency  com¬ 
ponents.  We  simulated  the  response  of  the  harbor  to  an  irregular  wave  train. 
Numerical  simulations  were  carried  out  for  an  incident  sea  state  characterized  by 
a  JONSWAP  spectrum  with  Hmo  =  3  m,  Tp  =  12  s  and  y  =  3.3.  Figure  41  shows  a 
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Figure  39.  Boussinesq  model  prediction  of  the  time-history  of  the  water-surface 
elevation  at  Gauge  5  for  a  natural  harbor  period  (7"  =  60  s) 


Figure  40.  Time-histories  of  total  and  long-period  (T  >  25  s)  component  of  water- 
surface  elevation  at  Gauges  1 ,  2,  3,  and  5  for  bichromatic  wave  train 
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Figure  41 .  3-D  view  of  irregular  wave  propagation  into  Barbers  Point  Harbor 

snapshot  of  waves  propagating  into  the  harbor.  The  predicted  wave  spectra  at 
Gauges  1  and  2  are  shown  in  Figure  42.  It  can  be  seen  that  nonlinear  interactions 
significantly  modified  the  wave  spectrum  at  Gauge  2{h  =  l  m),  with  the  cross- 
spectral  transfer  of  energy  to  both  lower  and  higher  harmonics.  The  wave  spectra 
for  the  four  gauges  located  inside  the  harbor  basin  are  plotted  in  Figure  43.  The 
short-period  wave  energy  is  much  smaller  inside  the  harbor  basin  because  the 
bathymetry  acts  to  refract  waves  away  from  the  harbor  entrance. 

The  wave  spectra  were  divided  into  short-period  ( T  <  25  s)  and  long-period 
components  ( T  >  25  s)  and  wave  heights  were  calculated  for  each  component. 
The  short-period  wave  heights  at  the  inside  gauges  varied  from  0.42  to  0.57  m, 
compared  to  3.1  m  at  outside  Gauge  2.  The  long-wave  energy  inside  the  harbor 
is,  however,  comparable  to  that  at  the  outside  gauge  with  heights  ranging  from 
0.38  to  0.46  m  for  the  inside  gauges,  compared  to  0.49  m  at  an  outside  Gauge  2. 
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Figure  42.  Predicted  wave  spectra  at  the  outside  Gauges  1  and  2  for  an  irregular 
sea  state  (Hmo  =  3  m,  Tp  =  12  s) 
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Figure  43.  Predicted  wave  spectra  at  gauges  inside  harbor  basin  (Gauges  3-6) 
for  an  irregular  sea  state  ( Hmo  =  3  m,  Tp  =  12  s) 
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Appendix  A 

Fourier  Series  Solutions  of 
Boussinesq  Equations 


As  the  height  of  surface  waves  in  intermediate  and  shallow-water  depths 
increases,  the  wave  profile  changes  from  a  sinusoidal  shape  to  an  asymmetric  one 
with  peaked  crests  and  broad  shallow  troughs.  The  Boussinesq  equations  are 
nonlinear  and  are  able  to  describe  the  change  in  wave  shape,  provided  the  wave 
height  and  period  are  within  nonlinear  and  dispersive  limits  of  the  equations. 
However,  it  is  important  that  the  velocity  and  flux  time-histories  imposed  along 
the  wave  generation  boundaries  of  the  numerical  model  be  consistent  with  the 
equations  that  are  solved  within  the  computational  domain.  In  BOUSS-2D,  the 
Fourier  approximation  method  of  Rienecker  and  Fenton  (1981)1  has  been  used  to 
derive  nonlinear  boundary  conditions  for  periodic  waves  in  water  of  constant 
depth.  The  one-dimensional  form  of  the  weakly  nonlinear  form  of  the  Boussinesq 
equations  (Equations  4  to  6)  for  water  of  constant  depth  can  be  written  as: 
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For  periodic  waves,  the  partial  differential  equations  can  be  transformed  into  a 
set  of  coupled  nonlinear  ordinary  differential  equations  in  terms  of  coordinate 
system,  £  =  x-C  t,  moving  at  the  phase  speed  of  the  waves,  C,  and  integrated 
once  to  yield: 
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1  References  cited  in  Appendices  A-E  are  listed  in  the  References  at  the  end  of  the  main  text. 
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Al 


where  Q  is  the  volume  flux  and  R  is  the  Bernoulli  constant.  The  velocity  ua  can  be 
expanded  as  a  Fourier  series: 


A2 


N 

Ua  =50  +2  JkBj  C0S(^)  (A5) 

j= 1 

where  N  is  the  number  of  Fourier  components,  Bj  are  the  Fourier  coefficients,  and 
k  is  the  wave  number.  To  solve  the  problem  numerically,  the  free  surface  eleva¬ 
tion  is  discretized  into  N+l  equally  spaced  points  over  half  a  wavelength,  i.e., 

T\n=*\fem),m=OX ,~.,N  (A6) 


where  E,m  =  mLUN  and  L  is  the  wavelength.  Equations  A3  to  A4  are  then  evalu¬ 
ated  at  points  over  half  a  wavelength  to  yield  a  system  of  nonlinear  algebraic 
equations: 
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The  above  2N+2  equations  involve  2N  +  5  unknowns  Ty(j  =  0, ...,  N),  Bj  (J  = 
0,  . . .,  N),  k,  Q,  R,  so  three  additional  equations  are  needed.  These  can  be  obtained 
from  the  wave  height,  H,  wave  period,  T,  the  mean  water  level  as: 
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N- 1 

R0  +1V  +  2IX  =  0 

(All) 

noting  that  C  =  -B0  for  a  zero-mean  Eulerian  velocity.  A  Newton-Raphson  pro¬ 
cedure  (see  Rienecker  and  Fenton  1981)  is  used  to  solve  the  system  of  equations 
(A7  to  A1 1)  for  the  unknown  values  of  the  free  surface  displacement  at  the 
collocation  points,  the  Fourier  coefficients,  the  wave  number,  the  phase  speed, 
and  the  constants  Q  and  R. 
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Appendix  B 

Description  of  Ocean  Wave 
Spectra 


A  number  of  parametric  shapes  have  been  proposed  to  describe  the  frequency 
distribution  of  wave  energy.  The  spectral  shapes  were  derived  from  long-term 
field  wave  measurements  and  depend  on  factors  such  as  wind  duration,  the 
distance  or  fetch  over  which  the  wind  blows,  and  water  depth. 


Pierson-Moskowitz  Spectrum 

For  fully-developed  sea  states  in  deep  water  where  there  is  a  local  balance 
between  momentum  transfer  from  the  wind  and  wave  breaking/nonlinear  cross- 
spectral  energy  transfer  processes,  Pierson-Moskowitz  (1964)1  proposed  the 
following  wave  spectrum: 


S(f)  = 


«g2 

(2n )4/5 


exp 


f  v* 
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where  a  =  0.0081  is  Phillips’  constant,  g  is  the  gravitational  acceleration,  = 
0.74,  fp-  g/2nU  i9.5,  and  U  19.5  is  the  wind  speed  at  19.5  m  above  the  mean  sea 
level. 


Bretschneider  Spectrum 

The  Bretschneider  spectrum  (Bretschneider  1959)  has  the  same  shape  as  the 
Pierson-Moskowitz  spectrum  but  is  defined  in  terms  of  the  significant  wave 
height,  Hs,  and  spectral  peak  frequency,^,,  instead  of  the  wind  speed.  It  can  be 
written  as: 


1 


References  cited  in  this  appendix  are  listed  in  the  References  at  the  end  of  the  main  text. 
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JONSWAP  Spectrum 

Based  on  an  extensive  analysis  of  data  from  the  Joint  North  Sea  Wave 
Observation  Project  (JONSWAP),  Hasselmann  et  al.  (1973)  proposed  the 
following  modified  form  of  the  Pierson-Moskowitz  spectrum  to  account  for  fetch- 
limited  conditions: 
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F  is  the  fetch  distance,  y  is  a  spectral  peak  enhancement  factor,  and  Uw  is  the 
wind  speed  at  10  m  above  the  mean  sea  level.  Figure  B1  shows  a  comparison  of  a 
JONSWAP  spectrum  with  y=  3.3  with  a  Bretschneider  spectrum  for  a  sea  state. 
The  JONSWAP  spectrum  reduces  to  the  Bretschneider  spectrum  when  y  is  equal 
to  1.0. 

TMA  Spectrum 

The  TMA  spectrum  is  a  modified  version  of  the  JONSWAP  spectrum  for 
water  of  finite  depth.  It  was  proposed  by  Bouws,  Gunther,  and  Vincent  (1985) 
based  on  Kitaigorodskii,  Krasitskii,  and  Zaslavskii  (1975)  frequency-dependent 
factor  for  the  equilibrium  range  of  the  wave  spectrum  in  shallow-water  and 
validated  with  data  from  three  field  studies  (Texel,  MARSEN,  and  ARSLOE).  It 
can  be  written  as: 
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Figure  B1 .  Comparison  of  Bretschneider  and  JONSWAP  (y  =  3.3)  spectra  for  a 
sea  state  with  Hmo  =1  m,  7p  =  10s 
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where  Mffi)  is  a  function  that  expresses  the  effect  of  finite  water  depth  and  is 
given  by: 
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The  frequency  factor  CD/,  =  2nf  ~J(h/  g)  and  R((£>h)  is  obtained  from  the  iterative 
solution  of  the  linear  dispersion  relation: 


i?(coA)tanh[(D^(coA)]  =  1 
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Ochi-Hubble  Spectrum 


The  Ochi-Hubble  spectrum  (Ochi  and  Hubble  1976)  is  a  six-parameter 
spectrum  proposed  for  bimodal  sea  states  consisting  of  swell  and  local  wind¬ 
generated  components.  It  can  be  written  as: 
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where  T  is  the  gamma  function.  Hs,fp,  and  /denote  the  significant  wave  height, 
spectral  peak  frequency,  and  spectral  shape  factor  respectively  for  the  two  sea 
state  components.  When  y  is  equal  to  1 ,  the  Ochi-Hubble  spectrum  reduces  to  the 
Bretschneider  spectrum. 
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Appendix  C 

Directional  Wave  Spreading 
Functions 


The  directional  spreading  function,  D(Q),  describes  the  directional  distribution 
of  wave  energy  in  irregular  multidirectional  sea  states.  It  can  be  quantified  in 
terms  of  the  principal  direction  of  wave  propagation  Qp,  and  the  directional  spread 
or  standard  deviation  of  the  spreading  function,  c@,  which  is  defined  as: 

n  r  0„+7t/2  ~ 

=  l%nD{Q)^P)dQ  (C1) 

A  number  of  parametric  shapes  have  been  proposed  to  describe  the  directional 
spreading  function  including  the  cosine-power,  the  circular  normal,  and  wrapped- 
normal  distributions.  These  are  described  in  the  following  paragraphs. 


Cosine-Power  Spreading  Function 

The  cosine-power  function  is  an  extended  version  of  the  cosine-squared  direc¬ 
tional  distribution  initially  proposed  by  St.  Denis  and  Pierson  (1953)1  and  can  be 
written  as: 

£>(Q)  =  cos^(e-e  )  for  10-0,  |  <  ti/2  (C2) 

V7ir(s  +  l/2)  p 

where  T  is  the  gamma  function.  The  parameter  s  is  an  index  describing  the  degree 
of  directional  spreading  with  s  — »  °°  representing  a  unidirectional  wave  field. 


Circular-Normal  Spreading  Function 

The  circular  normal  distribution  was  proposed  by  Borgman  (1969)  and  can  be 
written  as: 


References  cited  in  this  appendix  are  listed  in  the  References  at  the  end  of  the  main  text. 
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B(e)=w^jexp[,ICOS(0'0')] 

where  I0  is  the  modified  Bessel  function  of  the  first  kind  and  a  is  a  parameter 
describing  the  degree  of  directional  spreading  with  a  — >  °°  representing  a 
unidirectional  wave  field. 


Wrapped-Normal  Spreading  Function 

The  wrapped-normal  distribution  was  suggested  by  Mardia  (1972)  and  is 
given  by; 
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Figure  Cl  shows  a  plot  of  the  distributions  for  the  three  different  spreading  func¬ 
tion  formulations  corresponding  to  a  standard  deviation  Oe  of  25.5  deg.  The 
associated  spreading  indices  are  s  =  2  for  the  cosine-power  function  and  a  =  5.55 
for  the  circular-normal  distribution.  Thirty  components  (N=  30)  were  used  for  the 
wrapped-normal  distribution.  The  cosine-normal  and  wrapped-normal  distribu¬ 
tions  are  slightly  narrower  than  the  cosine-power  function  although  the  differences 
can  be  considered  to  be  minimal. 


0  (deg) 


Figure  Cl.  Comparison  of  the  cosine-power,  circular-normal  and  wrapped-normal 
distributions  with  a  standard  deviation  of  25.5  deg 
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Appendix  D 

B0USS-2D  File  Formats 


BOUSS-2D  currently  supports  two  file  formats.  The  grid  format  (.grd)  is  used 
to  define  the  spatial  variation  of  scalar  or  vector  quantities  over  a  rectangular  grid. 
The  time  series  format  (.tsl)  is  used  for  one-dimensional  time-histories  of  scalar 
quantities.  The  files  consist  of  a  header  section  and  a  data  section.  The  header 
section  includes  both  mandatory  parameters  required  to  define  the  data  structure 
and  optional  parameters  that  provide  additional  information  on  the  data.  Comment 
lines  in  the  header  section  begin  with  the  pound  (#)  sign,  and  lines  with  parameter 
names  and  values  begin  with  the  colon  (:)  sign. 


Grid  File  Format  (.grd) 

The  grid  file  format  is  used  to  store  multiple  2-D  arrays  hj{x,y),  h2{x,y),  etc. 
The  variables  are  defined  over  a  rectangular  grid  with  the  origin  located  at  (xongm, 
y origin)  and  uniform  grid  spacings  Ax  and  Ay  in  the  x  and  y  directions  respectively. 
The  number  of  grid  points  are  Nx  and  Ny  in  the  x  and  y  directions  respectively. 
The  2-D  array  h(ij)  represents  the  value  of  h(xy)  at  x  =  xorigm  +  (z'-l)Ax  andy  = 
yorigin  +  (/-l)Ay.  The  mandatory  header  parameters  are: 

•  NGrid  X  =  number  of  grid  points  in  the  x  direction  (Nx) 

•  NGrid  Y  =  number  of  grid  points  in  the  y  direction  (Ny) 

•  Delta  X  =  grid  spacing  in  the  x  direction  (Ax) 

•  Delta  Y  =  grid  spacing  in  the  y  direction  (Ay) 

•  X  Origin  =  x-coordinate  of  grid  origin 

•  Y  Origin  =  y-coordinate  of  grid  origin 

•  N  Arrays  =  number  of  data  arrays  (1  for  scalar  quantities,  2  for  vector 
quantities) 

The  header  section  must  terminate  with  the  EndHeader  keyword.  Two- 
dimensional  data  are  then  written  out  after  the  header  section  as  continuous  data 
streams  in  row  first  order,  i.e.,  (( h(ij ),  /  =  1,  Nx),j  -  1,  Ny)  for  scalar  quantities  or 
i  =  1,  Nx),j  =  1 ,  Ny)  and  ((v(y),  i  =  1,  Nx),j  =  1,  Ny)  for  vector  quantities. 
A  sample  grid  file  is  provided  as  follows: 
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#  B0USS2D  Grid  File 

:FileName  bathy.grd 

:WrittenBy  Bouss2D  Version  1.0 

: CreationDate  Mon  Sep  18  12:40:25  2000 

# 


N_Arrays  1 

DataDescription (1)  Seabed  Elevation 
DataUnits (1)  meters 


# 


:NGrid_X 

4 

:NGrid  Y 

4 

:  Delta  X 

10.0 

:  Delta_Y 

10 . 0 

: X  Origin 

0.0 

:Y  Origin 
# 

0.0 

:  EndHeader 

15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8  15.8 


Time  Series  File  Format  (.tsl) 

The  time  series  file  format  is  used  to  store  time-histories  of  scalar  quantities 
u(t )  at  discrete  time  intervals  t,  =  fi  +  (/- 1 )  A/ ,  where  t\  is  the  start  time  and  At  is 
the  time-step.  The  number  of  data  points  in  the  time  record  is  N,.  Multiple  time- 
histories  (e.g.,  data  output  at  different  grid  locations)  can  be  stored  in  the  same 
file.  The  mandatory  header  parameters  are: 

•  NTimeSteps  =  number  of  time-steps  (N,) 

•  TimeStep  =  time-step  (At) 

•  StartTime  =  start  time  of  record  (s) 

•  NDatasets  =  number  of  time  series  records 

The  header  section  terminates  with  the  EndHeader  keyword.  Data  are  then 
written  out  after  the  header  section  as  a  single  column  vector  for  one  time  record 
or  multiple  columns  for  multiple  time  records.  A  sample  time  series  file  is 
provided  as  follows: 

•  BOUSS2D  Time  Series  File 

:FileName  test. tsl 

:WrittenBy  Bouss2D  Version  1.0 


: CreationDate  Mon 

.u 

Sep  18  12: 

:  4 0 : 25  2000 

: NDataSets 

4 

: DataDescription (1) 

Surface 

Elevation 

: DataUnits (1) 

meters 

: DataDescription (2 ) 

Surface 

Elevation 

: DataUnits (2 ) 

meters 

: DataDescription (3 ) 

Surface 

Elevation 

: DataUnits (3 ) 

: DataDescription {4 ) 

: DataUnits (4 ) 

# 

: StartTime 

meters 

Surface 

meters 

Elevation 

0.0000 

: TimeStep 

0.1500 

: NTimeSteps 

5334 

# 
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:x_grid (1) 

270.0000 

:y_grid (1) 

1030.0000 

:x_grid (2) 

500.0000 

:y_grid (2) 

1030.0000 

:x_grid (3 ) 

820.0000 

:y_grid (3) 

1030.0000 

:x  grid(4) 

1000.0000 

:y__grid  (4 ) 

1030.0000 

it 

:  EndHeader 

0 . OOOOOOE+OO 

-1 . 772061E- 15 

-1 . 195287E-13 

-2 . 679531E-13 

- 1 . 777521E- 12 

-3 . 01248 8E- 12 

-1 . 2243 17E- 11 

- 1 . 867146E- 11 

-5 . 994772E- 11 

-8 . 602132E-11 

-2 . 376610E- 10 

-3 . 269985E-10 

-  8 . 099900E- 10 

-1 . 079451E-09 

-2 . 454132E- 09 

-3 . 187200E- 09 

-6 . 7496  91E- 09 

- 8 . 576622E- 09 

-1 . 234493E-14  -4 . 485098E- 14 
-5 . 388106E-13  -1 . 005685E- 12 
-4 . 936255E-12  -7 . 8662 51E- 12 
-2 . 797059E-11  -4 . 124122E- 11 
-1 .219879E-10  - 1 . 711262E- 10 
-4 . 460105E- 10  -6. 033756E-10 
-1 . 428627E-09  - 1 . 87833 8E- 09 
-4 . 115422E- 09  -5 . 284574E- 09 
-1. 084376E-08  -1 . 364400E- 08 
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Appendix  E 
Utility  Programs 


GEN_DAMP 

GENDAMP  generates  a  damping  file  given  the  location  of  the  damping 
layers  (north,  south,  east,  or  west  boundaries  of  the  grid),  the  width  of  the 
damping  layer,  and  the  nondimensional  damping  strength  at  the  end  of  the 
damping  layer,  /v  The  damping  values  are  initially  set  to  zero  over  the  entire 
grid.  The  program  then  calculates  the  spatial  variation  of  the  damping  values  over 
the  width  of  the  damping  layers  using  a  quadratic  function,  e.g., 


v«d(x>y) 


x-xn 


K 


where  x0  and  xw  are  the  coordinates  of  the  beginning  and  end  of  the  damping 
layer.  The  program  is  run  interactively  from  an  MS-DOS  command  prompt 
window.  An  example  is  as  follows: 

C : \bouss2d\bin\gen_damp 

bathy  #  name  of  input  bathymetry  file  (.grd) 
damp  #  name  of  output  damping  file  (.grd) 

0  #  width  of  damping  layer  -  North  (m)  [0.0] 

0  #  non-dimensional  damping  strength  (0-1)  [0.0] 

100  #  width  of  damping  layer  -  East  (m)  [0.0] 

1  #  non-dimensional  damping  strength  (0-1)  [0.0] 

0  #  width  of  damping  layer  -  South  (m)  [0.0] 

0  #  non-dimensional  damping  strength  (0-1)  [0.0] 

0  #  width  of  damping  layer  -  West  (m)  [0.0] 

0  #  non-dimensional  damping  strength  (0-1)  [0.0] 

To  effectively  damp  out  waves  at  open  boundaries,  a  damping  layer  half  a  wave¬ 
length  wide  with  a  damping  strength  of  1 .0  should  be  used.  The  damping  file  is 
written  out  as  an  ASCII  file  in  the  grid  file  format  described  in  Appendix  D. 


MAP  POROSITY 


MAPPOROSITY  creates  a  porosity  grid  file  given  the  boundaries  of  porous 
regions  within  the  computational  domain  as  a  set  of  discrete  (x,y)  points.  The 
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porosity  values  are  initially  set  to  1.0  over  the  entire  grid.  The  program  then  uses 
the  “Winding-Number”  algorithm  of  Gordkin  and  Pulli  (1984)1  to  set  the  porosity 
values  within  the  porous  region  boundaries  to  values  specified  by  the  user.  The 
(x,y)  coordinates  of  the  porous  region  boundaries  are  stored  as  two  columns  of 
data  in  an  ASCII  file.  The  program  is  run  interactively  from  an  MS-DOS 
command  prompt  window.  An  example  is  as  follows: 


C:  \bouss2d\bin\map_porosity 

bathy  #  name  of  input  bathymetry  file  (.grd) 

porosity  #  name  of  output  porosity  file  (.grd) 

2  #  number  of  porous  regions  [1] 

bwl.xy  #  name  of  boundary  file  for  porous  region  #1 

0.4  #  porosity  (0-1)  [0.4] 

bw2.xy  #  name  of  boundary  file  for  porous  region  #2 

0.4  #  porosity  (0-1)  [0.4] 


[.xy] 

[.xy] 


The  porosity  file  is  written  out  as  an  ASCII  file  in  the  grid  file  format  described  in 
Appendix  D. 


1  References  cited  in  this  appendix  are  listed  in  the  References  at  the  end  of  the  main  text. 
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